#ifndef GITHASH_PP
#define GITHASH_PP "unknown"
#endif
#ifndef GITDATE_PP
#define GITDATE_PP "unknown"
#endif
#ifndef BUILDDATE_PP
#define BUILDDATE_PP "unknown"
#endif
MODULE COMP_FUNCTIONS

! I/O + OS functions

USE PRECISION_PARAMETERS
IMPLICIT NONE

CONTAINS

REAL(EB) FUNCTION CURRENT_TIME()  ! Returns the wall clock time in seconds.
USE MPI
CURRENT_TIME = MPI_WTIME()
END FUNCTION CURRENT_TIME


SUBROUTINE GET_DATE(DATE)

INTEGER :: DATE_TIME(8)
CHARACTER(10) :: BIG_BEN(3),MONTH
CHARACTER(LABEL_LENGTH), INTENT(OUT) :: DATE

CALL DATE_AND_TIME(BIG_BEN(1),BIG_BEN(2),BIG_BEN(3),DATE_TIME)

SELECT CASE(DATE_TIME(2))
   CASE(1)
      MONTH='January'
   CASE(2)
      MONTH='February'
   CASE(3)
      MONTH='March'
   CASE(4)
      MONTH='April'
   CASE(5)
      MONTH='May'
   CASE(6)
      MONTH='June'
   CASE(7)
      MONTH='July'
   CASE(8)
      MONTH='August'
   CASE(9)
      MONTH='September'
   CASE(10)
      MONTH='October'
   CASE(11)
      MONTH='November'
   CASE(12)
      MONTH='December'
END SELECT

WRITE(DATE,'(A,I3,A,I4,2X,I2.2,A,I2.2,A,I2.2)') TRIM(MONTH),DATE_TIME(3),', ',DATE_TIME(1), &
                                                DATE_TIME(5),':',DATE_TIME(6),':',DATE_TIME(7)
END SUBROUTINE GET_DATE


SUBROUTINE SHUTDOWN(MESSAGE,PROCESS_0_ONLY)

! Stops the code gracefully after writing a message

USE GLOBAL_CONSTANTS, ONLY: MYID,LU_ERR,CHID,SETUP_STOP,STOP_STATUS,N_MPI_PROCESSES
CHARACTER(*) MESSAGE
LOGICAL, INTENT(IN), OPTIONAL :: PROCESS_0_ONLY
CHARACTER(20) :: FMT

IF (PRESENT(PROCESS_0_ONLY)) THEN
   WRITE(FMT,'("(/A,A,I",I0,",A,A,A)")') CEILING(LOG10(REAL(N_MPI_PROCESSES+1,EB))) + 1
   IF (.NOT.PROCESS_0_ONLY) WRITE(LU_ERR,FMT) TRIM(MESSAGE),' (MPI Process:',MYID,', CHID: ',TRIM(CHID),')'
   IF (PROCESS_0_ONLY .AND. MYID==0) WRITE(LU_ERR,'(/A,A,A,A)') TRIM(MESSAGE),' (CHID: ',TRIM(CHID),')'
ELSE
   IF (MYID==0) WRITE(LU_ERR,'(/A,A,A,A)') TRIM(MESSAGE),' (CHID: ',TRIM(CHID),')'
ENDIF

STOP_STATUS = SETUP_STOP

END SUBROUTINE SHUTDOWN


SUBROUTINE SYSTEM_MEM_USAGE(VALUE_RSS)

! This routine returns the memory used by the given process at the time of the call. It only works on a Linux machine because
! it searches for the system file called '/proc/PID/status' and reads the VmRSS value (kB).
! The DEVC output QUANTITY 'RAM' outputs this value in units of MB.
! Normally, this routine just returns 0 because it is non-standard Fortran and we have disconnected the GETPID call.
! To invoke it, uncomment the GETPID statement, and the 'USE IFPORT' statement if compiling with Intel Fortran.

!USE IFPORT  ! Intel Fortran extension library. This is needed for GETPID.
INTEGER, INTENT(OUT) :: VALUE_RSS
CHARACTER(200):: FILENAME=' '
CHARACTER(80) :: LINE
CHARACTER(8)  :: PID_CHAR=' '
INTEGER :: PID
LOGICAL :: IFXST

VALUE_RSS=0  ! return zero if the memory file is not found

PID = 0
!PID = GETPID() ! GETPID is non-standard Fortran. Leave it commented out for the moment.

WRITE(PID_CHAR,'(I8)') PID
FILENAME = '/proc/'//TRIM(ADJUSTL(PID_CHAR))//'/status'

INQUIRE(FILE=FILENAME,EXIST=IFXST)
IF (.NOT.IFXST) then
   WRITE (*,*) 'WARNING: Memory system file does not exist'
  RETURN
ENDIF

OPEN(1000,FILE=FILENAME,ACTION='READ')
DO
   READ (1000,'(A)',END=120) LINE
   IF (LINE(1:6)=='VmRSS:') THEN
      READ (LINE(7:),*) VALUE_RSS
      EXIT
   ENDIF
ENDDO
120 CONTINUE
CLOSE(1000)

END SUBROUTINE SYSTEM_MEM_USAGE


SUBROUTINE GET_INPUT_FILE ! Read the argument after the command
USE GLOBAL_CONSTANTS, ONLY: FN_INPUT
IF (FN_INPUT=='null') CALL GET_COMMAND_ARGUMENT(1,FN_INPUT)
END SUBROUTINE GET_INPUT_FILE


INTEGER FUNCTION GET_FILE_NUMBER()
USE GLOBAL_CONSTANTS, ONLY: MYID,FILE_COUNTER
FILE_COUNTER(MYID) = FILE_COUNTER(MYID) + 1
GET_FILE_NUMBER = FILE_COUNTER(MYID)
END FUNCTION GET_FILE_NUMBER


SUBROUTINE CHECKREAD(NAME,LU,IOS)

! Look for the namelist variable NAME and then stop at that line.

USE GLOBAL_CONSTANTS, ONLY: INPUT_FILE_LINE_NUMBER
INTEGER :: II
INTEGER, INTENT(OUT) :: IOS
INTEGER, INTENT(IN) :: LU
CHARACTER(4), INTENT(IN) :: NAME
CHARACTER(80) TEXT
IOS = 1

READLOOP: DO
   READ(LU,'(A)',END=10) TEXT
   INPUT_FILE_LINE_NUMBER = INPUT_FILE_LINE_NUMBER + 1
   TLOOP: DO II=1,72
      IF (TEXT(II:II)/='&' .AND. TEXT(II:II)/=' ') EXIT TLOOP
      IF (TEXT(II:II)=='&') THEN
         IF (TEXT(II+1:II+4)==NAME) THEN
            BACKSPACE(LU)
            IOS = 0
            EXIT READLOOP
         ELSE
            CYCLE READLOOP
         ENDIF
      ENDIF
   ENDDO TLOOP
ENDDO READLOOP

10 RETURN
END SUBROUTINE CHECKREAD


SUBROUTINE SCAN_INPUT_FILE(LU,IOS,TEXT)

! Look for odd or illegal characters in the input file.

INTEGER, INTENT(OUT) :: IOS
INTEGER, INTENT(IN) :: LU
CHARACTER(80), INTENT(OUT) :: TEXT

IOS = 1

READLOOP: DO
   READ(LU,'(A)',END=10) TEXT
   IF (IACHAR(TEXT(1:1))==13 .AND. TEXT(2:2)=='&') THEN
      IOS = 0
      EXIT READLOOP
   ENDIF
ENDDO READLOOP

10 RETURN
END SUBROUTINE SCAN_INPUT_FILE


SUBROUTINE APPEND_FILE(LU,N_TEXT_LINES,T)

! Read through a file and stop when T exceeds the last time read

USE GLOBAL_CONSTANTS, ONLY: CLIP_RESTART_FILES
INTEGER, INTENT(IN) :: LU,N_TEXT_LINES
INTEGER :: ITER
REAL(EB) :: FILE_TIME
REAL(EB), INTENT(IN) :: T

READ_TEXT_LINES: DO ITER=1,N_TEXT_LINES
   READ(LU,*)
ENDDO READ_TEXT_LINES

READ_LOOP: DO
   READ(LU,*,END=10) FILE_TIME
   IF (FILE_TIME>T .AND. CLIP_RESTART_FILES) THEN
      BACKSPACE(LU)
      EXIT READ_LOOP
   ENDIF
ENDDO READ_LOOP

RETURN
10 BACKSPACE(LU)

END SUBROUTINE APPEND_FILE


SUBROUTINE CHECK_XB(XB)
! Reorder an input sextuple XB if needed
REAL(EB) :: DUMMY,XB(6)
INTEGER  :: I
DO I=1,5,2
   IF (XB(I)>XB(I+1)) THEN
      DUMMY   = XB(I)
      XB(I)   = XB(I+1)
      XB(I+1) = DUMMY
   ENDIF
ENDDO
END SUBROUTINE CHECK_XB


SUBROUTINE CHANGE_UNITS(QUANTITY,UNITS,SPATIAL_STATISTIC,TEMPORAL_STATISTIC,MYID,LU_ERR)

! Change the units of the DEVC output if it is an integrated quantity

CHARACTER(LABEL_LENGTH), INTENT(IN) :: QUANTITY,SPATIAL_STATISTIC,TEMPORAL_STATISTIC
INTEGER, INTENT(IN) :: MYID,LU_ERR
CHARACTER(LABEL_LENGTH), INTENT(INOUT) ::UNITS
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
INTEGER :: I,UNIT_L

UNIT_L = LEN(TRIM(UNITS))

I = 1
SELECT CASE (SPATIAL_STATISTIC)
   CASE('VOLUME INTEGRAL')
      I = INDEX(UNITS,'/m3')
      IF (I/=0) WRITE(UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+3:UNIT_L)
      IF (TRIM(UNITS)=='1/s') THEN
         UNITS = 'm3/s'
         I = 1
      ENDIF
      IF (TRIM(UNITS)=='Pa') THEN
         UNITS = 'J'
         I = 1
      ENDIF
   CASE('MASS INTEGRAL')
      I = INDEX(UNITS,'/kg')
      IF (I/=0) WRITE(UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+3:UNIT_L)
   CASE('AREA INTEGRAL','SURFACE INTEGRAL')
      I = INDEX(UNITS,'/m2')
      IF (I/=0) WRITE(UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+3:UNIT_L)
   CASE ('VOLUME')
      UNITS = 'm3'
      I=1
   CASE ('AREA','SURFACE AREA')
      UNITS = 'm2'
      I=1
END SELECT

IF (I==0) THEN
   WRITE(MESSAGE,'(A,A,A,A)')  'WARNING: Problem with units compatibility of SPATIAL_STATISITIC ',TRIM(SPATIAL_STATISTIC), &
                               ' with the QUANTITY ',TRIM(QUANTITY)
   IF (MYID==0) WRITE(LU_ERR,'(A)') TRIM(MESSAGE)
ENDIF

I = 1
SELECT CASE (TEMPORAL_STATISTIC)
   CASE('TIME INTEGRAL')
      I = INDEX(UNITS,'/s')
      IF (I/=0) THEN
         WRITE(UNITS,'(A,A)') UNITS(1:I-1),UNITS(I+2:UNIT_L)
      ELSE
         I = INDEX(UNITS,'kW')
         IF (I==1) WRITE(UNITS,'(A,A)') 'kJ',UNITS(3:UNIT_L)
         IF (I>1)  WRITE(UNITS,'(A,A,A)') UNITS(1:I-1),'kJ',UNITS(I+2:UNIT_L)
      ENDIF
      IF (QUANTITY=='FED') I=1
END SELECT

IF (I==0) THEN
   WRITE(MESSAGE,'(A,A,A,A)')  'WARNING: Problem with units compatibility of TEMPORAL_STATISTIC ',TRIM(TEMPORAL_STATISTIC), &
                               ' with the QUANTITY ',TRIM(QUANTITY)
   IF (MYID==0) WRITE(LU_ERR,'(A)') TRIM(MESSAGE)
ENDIF

END SUBROUTINE CHANGE_UNITS


END MODULE COMP_FUNCTIONS



MODULE MEMORY_FUNCTIONS

USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE COMP_FUNCTIONS, ONLY : SHUTDOWN
IMPLICIT NONE

CONTAINS


SUBROUTINE ChkMemErr(CodeSect,VarName,IZERO)

! Memory checking routine

CHARACTER(*), INTENT(IN) :: CodeSect, VarName
INTEGER IZERO
CHARACTER(MESSAGE_LENGTH) MESSAGE

IF (IZERO==0) RETURN

WRITE(MESSAGE,'(4A)') 'ERROR: Memory allocation failed for ', TRIM(VarName),' in the routine ',TRIM(CodeSect)
CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)

END SUBROUTINE ChkMemErr



FUNCTION REALLOCATE(P,N1,N2)

! Resize the array P

REAL(EB), POINTER, DIMENSION(:) :: P, REALLOCATE
INTEGER, INTENT(IN) :: N1,N2
INTEGER :: NOLD, IERR
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
ALLOCATE(REALLOCATE(N1:N2), STAT=IERR)
IF (IERR /= 0) THEN
   WRITE(MESSAGE,'(A)') 'ERROR: Memory allocation failed in REALLOCATE'
   CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
ENDIF
IF (.NOT. ASSOCIATED(P)) RETURN
NOLD = MIN(SIZE(P), N2-N1+1)
REALLOCATE(N1:NOLD+N1-1) = P(N1:NOLD+N1-1)  ! Restore the contents of the reallocated array
DEALLOCATE(P)

END FUNCTION REALLOCATE


FUNCTION REALLOCATE_CHARACTER_ARRAY(P,CLEN,N1,N2)

! Resize the character array P

INTEGER, INTENT(IN) :: N1,N2,CLEN
CHARACTER(CLEN), POINTER, DIMENSION(:) :: P, REALLOCATE_CHARACTER_ARRAY
INTEGER :: NOLD, IERR
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
ALLOCATE(REALLOCATE_CHARACTER_ARRAY(N1:N2), STAT=IERR)
IF (IERR /= 0) THEN
   WRITE(MESSAGE,'(A)') 'ERROR: Memory allocation failed in REALLOCATE_CHARACTER_ARRAY'
   CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
ENDIF
IF (.NOT. ASSOCIATED(P)) RETURN
NOLD = MIN(SIZE(P), N2-N1+1)
REALLOCATE_CHARACTER_ARRAY(N1:NOLD+N1-1) = P(N1:NOLD+N1-1)  ! Restore the contents of the reallocated array
DEALLOCATE(P)

END FUNCTION REALLOCATE_CHARACTER_ARRAY


FUNCTION REALLOCATE2D(P,M1,M2,N1,N2)

! Resize the 2D array P

REAL(EB), POINTER, DIMENSION(:,:) :: P, REALLOCATE2D
INTEGER, INTENT(IN) :: M1,M2,N1,N2
INTEGER :: MOLD,NOLD,IERR
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
ALLOCATE(REALLOCATE2D(M1:M2,N1:N2), STAT=IERR)
IF (IERR /= 0) THEN
   WRITE(MESSAGE,'(A)') 'ERROR: Memory allocation failed in REALLOCATE2D'
   CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
ENDIF
IF (.NOT. ASSOCIATED(P)) RETURN
MOLD = MIN(SIZE(P,DIM=1), M2-M1+1)
NOLD = MIN(SIZE(P,DIM=2), N2-N1+1)
REALLOCATE2D(M1:MOLD+M1-1,N1:NOLD+N1-1) = P(M1:MOLD+M1-1,N1:NOLD+N1-1)  ! Restore the contents of the reallocated array
DEALLOCATE(P)

END FUNCTION REALLOCATE2D


SUBROUTINE RE_ALLOCATE_STRINGS(NM)

CHARACTER(MESH_STRING_LENGTH), ALLOCATABLE, DIMENSION(:) :: DUMMY
INTEGER, INTENT(IN) :: NM
TYPE (MESH_TYPE), POINTER :: M=>NULL()

M=>MESHES(NM)
ALLOCATE(DUMMY(1:M%N_STRINGS))
DUMMY = M%STRING
DEALLOCATE(M%STRING)
ALLOCATE(M%STRING(M%N_STRINGS_MAX+100))
M%STRING(1:M%N_STRINGS) = DUMMY(1:M%N_STRINGS)
M%N_STRINGS_MAX = M%N_STRINGS_MAX+100
DEALLOCATE(DUMMY)

END SUBROUTINE RE_ALLOCATE_STRINGS



SUBROUTINE COMPUTE_ONE_D_STORAGE_DIMENSIONS(SURF_INDEX)

USE GLOBAL_CONSTANTS, ONLY: N_TRACKED_SPECIES
INTEGER, INTENT(IN) :: SURF_INDEX
TYPE (SURFACE_TYPE), POINTER :: SF

SF => SURFACE(SURF_INDEX)

SF%ONE_D_REALS_ARRAY_SIZE( 1) = N_TRACKED_SPECIES            ! MASSFLUX_SPEC
SF%ONE_D_REALS_ARRAY_SIZE( 2) = SF%N_MATL                    ! MASSFLUX_MATL
SF%ONE_D_REALS_ARRAY_SIZE( 3) = N_TRACKED_SPECIES            ! MASSFLUX
SF%ONE_D_REALS_ARRAY_SIZE( 4) = SF%N_CELLS_MAX+1             ! X
SF%ONE_D_REALS_ARRAY_SIZE( 5) = SF%N_CELLS_MAX+2             ! TMP
SF%ONE_D_REALS_ARRAY_SIZE( 6) = (SF%N_CELLS_MAX+2)*SF%N_MATL ! RHO
SF%ONE_D_REALS_ARRAY_SIZE( 7) = SF%N_LAYERS                  ! LAYER_THICKNESS
SF%ONE_D_REALS_ARRAY_SIZE( 8) = SF%NRA*SF%NSB                ! ILW
SF%ONE_D_REALS_ARRAY_SIZE( 9) = SF%NSB                       ! IL
SF%ONE_D_REALS_ARRAY_SIZE(10) = N_TRACKED_SPECIES            ! ZZ_G
SF%ONE_D_REALS_ARRAY_SIZE(11) = N_TRACKED_SPECIES            ! ZZ_F
SF%ONE_D_REALS_ARRAY_SIZE(12) = N_TRACKED_SPECIES            ! RHO_D_F
SF%ONE_D_REALS_ARRAY_SIZE(13) = N_TRACKED_SPECIES            ! RHO_D_DZDN_F

SF%ONE_D_INTEGERS_ARRAY_SIZE(1) = SF%N_LAYERS           ! N_LAYER_CELLS

SF%N_ONE_D_STORAGE_REALS    = N_ONE_D_SCALAR_REALS    + SUM(SF%ONE_D_REALS_ARRAY_SIZE(1:13))
SF%N_ONE_D_STORAGE_INTEGERS = N_ONE_D_SCALAR_INTEGERS + SUM(SF%ONE_D_INTEGERS_ARRAY_SIZE(1:1))
SF%N_ONE_D_STORAGE_LOGICALS = N_ONE_D_SCALAR_LOGICALS

END SUBROUTINE COMPUTE_ONE_D_STORAGE_DIMENSIONS


SUBROUTINE COMPUTE_WALL_STORAGE_DIMENSIONS(SURF_INDEX)

USE GLOBAL_CONSTANTS, ONLY: N_LP_ARRAY_INDICES,N_SURFACE_DENSITY_SPECIES
INTEGER, INTENT(IN) :: SURF_INDEX
TYPE (SURFACE_TYPE), POINTER :: SF

SF => SURFACE(SURF_INDEX)

SF%WALL_REALS_ARRAY_SIZE(1) = N_LP_ARRAY_INDICES            ! A_LP_MPUA
SF%WALL_REALS_ARRAY_SIZE(2) = N_LP_ARRAY_INDICES            ! LP_CPUA
SF%WALL_REALS_ARRAY_SIZE(3) = N_LP_ARRAY_INDICES            ! LP_MPUA
SF%WALL_REALS_ARRAY_SIZE(4) = N_SURFACE_DENSITY_SPECIES     ! AWM_AEROSOL
SF%WALL_REALS_ARRAY_SIZE(5) = SF%NVEG_L                     ! VEG_FUELMASS_L
SF%WALL_REALS_ARRAY_SIZE(6) = SF%NVEG_L                     ! VEG_MOISTMASS_L
SF%WALL_REALS_ARRAY_SIZE(7) = SF%NVEG_L+1                   ! VEG_TMP_L

SF%N_WALL_STORAGE_REALS    = N_WALL_SCALAR_REALS    + SUM(SF%WALL_REALS_ARRAY_SIZE(1:7))    + SF%N_ONE_D_STORAGE_REALS
SF%N_WALL_STORAGE_INTEGERS = N_WALL_SCALAR_INTEGERS                                         + SF%N_ONE_D_STORAGE_INTEGERS
SF%N_WALL_STORAGE_LOGICALS = N_WALL_SCALAR_LOGICALS                                         + SF%N_ONE_D_STORAGE_LOGICALS

END SUBROUTINE COMPUTE_WALL_STORAGE_DIMENSIONS


SUBROUTINE COMPUTE_CFACE_STORAGE_DIMENSIONS(SURF_INDEX)

USE GLOBAL_CONSTANTS, ONLY: N_SURFACE_DENSITY_SPECIES !,N_LP_ARRAY_INDICES
INTEGER, INTENT(IN) :: SURF_INDEX
TYPE (SURFACE_TYPE), POINTER :: SF

SF => SURFACE(SURF_INDEX)

!SF%CFACE_REALS_ARRAY_SIZE(1) = N_LP_ARRAY_INDICES            ! A_LP_MPUA
!SF%CFACE_REALS_ARRAY_SIZE(2) = N_LP_ARRAY_INDICES            ! LP_CPUA
!SF%CFACE_REALS_ARRAY_SIZE(3) = N_LP_ARRAY_INDICES            ! LP_MPUA
SF%CFACE_REALS_ARRAY_SIZE(4) = N_SURFACE_DENSITY_SPECIES     ! AWM_AEROSOL
!SF%CFACE_REALS_ARRAY_SIZE(5) = SF%NVEG_L                     ! VEG_FUELMASS_L
!SF%CFACE_REALS_ARRAY_SIZE(6) = SF%NVEG_L                     ! VEG_MOISTMASS_L
!SF%CFACE_REALS_ARRAY_SIZE(7) = SF%NVEG_L+1                   ! VEG_TMP_L

SF%N_CFACE_STORAGE_REALS    = N_CFACE_SCALAR_REALS   + SUM(SF%CFACE_REALS_ARRAY_SIZE(1:7)) + SF%N_ONE_D_STORAGE_REALS
SF%N_CFACE_STORAGE_INTEGERS = N_CFACE_SCALAR_INTEGERS                                      + SF%N_ONE_D_STORAGE_INTEGERS
SF%N_CFACE_STORAGE_LOGICALS = N_CFACE_SCALAR_LOGICALS                                      + SF%N_ONE_D_STORAGE_LOGICALS

END SUBROUTINE COMPUTE_CFACE_STORAGE_DIMENSIONS


SUBROUTINE COMPUTE_PARTICLE_STORAGE_DIMENSIONS(LPC_INDEX)

INTEGER, INTENT(IN) :: LPC_INDEX
TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
TYPE (SURFACE_TYPE), POINTER :: SF

LPC => LAGRANGIAN_PARTICLE_CLASS(LPC_INDEX)
SF  => SURFACE(LPC%SURF_INDEX)

LPC%N_STORAGE_REALS    = N_PARTICLE_SCALAR_REALS    + SF%N_ONE_D_STORAGE_REALS
LPC%N_STORAGE_INTEGERS = N_PARTICLE_SCALAR_INTEGERS + SF%N_ONE_D_STORAGE_INTEGERS
LPC%N_STORAGE_LOGICALS = N_PARTICLE_SCALAR_LOGICALS + SF%N_ONE_D_STORAGE_LOGICALS

END SUBROUTINE COMPUTE_PARTICLE_STORAGE_DIMENSIONS


SUBROUTINE ALLOCATE_STORAGE(NM,SURF_INDEX,LPC_INDEX,WALL_INDEX,CFACE_INDEX,LP_INDEX,TAG,NEW_TAG)

! This routine allocates space in a storage array for either a lagrangian particle or a wall cell

USE GLOBAL_CONSTANTS, ONLY: SIGMA,TMPA4,RPI,T_BEGIN,DIRICHLET,INITIALIZATION_PHASE,PERIODIC_TEST,TMPA
LOGICAL, INTENT(IN), OPTIONAL :: NEW_TAG
INTEGER, INTENT(IN):: NM,SURF_INDEX
INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,CFACE_INDEX,LP_INDEX,LPC_INDEX,TAG
INTEGER :: N_NEW_STORAGE_SLOTS,STORAGE_INDEX,I1,I2,I,OLD_INDEX,STORAGE_INDEX_OLD,STORAGE_INDEX_OPEN,N_STORAGE_SLOTS_OLD,NAS
LOGICAL :: NEW_SLOT
TYPE (STORAGE_TYPE), POINTER :: OS
TYPE(SURFACE_TYPE), POINTER :: SF=>NULL()
TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
TYPE (LAGRANGIAN_PARTICLE_TYPE), POINTER :: LP
TYPE (LAGRANGIAN_PARTICLE_TYPE), ALLOCATABLE, DIMENSION(:) :: LP_DUMMY
TYPE (WALL_TYPE), ALLOCATABLE, DIMENSION(:) :: WALL_DUMMY
TYPE (CFACE_TYPE), ALLOCATABLE, DIMENSION(:) :: CFACE_DUMMY
TYPE (MESH_TYPE), POINTER :: M

M  => MESHES(NM)
SF => SURFACE(SURF_INDEX)

IF (PERIODIC_TEST==7) THEN
   N_NEW_STORAGE_SLOTS=2*(M%IBAR*M%JBAR + M%IBAR*M%KBAR + M%JBAR*M%KBAR)
ELSE
   N_NEW_STORAGE_SLOTS=1000
ENDIF

IF (PRESENT(WALL_INDEX)) THEN

   OS => MESHES(NM)%WALL_STORAGE(SURF_INDEX)

   ! Check if there are enough wall cells

   IF (WALL_INDEX>M%N_WALL_CELLS_DIM) THEN
      ALLOCATE(WALL_DUMMY(0:M%N_WALL_CELLS_DIM))
      WALL_DUMMY = M%WALL
      DEALLOCATE(M%WALL)
      ALLOCATE(M%WALL(0:M%N_WALL_CELLS_DIM+N_NEW_STORAGE_SLOTS))
      M%WALL(0:MESHES(NM)%N_WALL_CELLS_DIM) = WALL_DUMMY
      DEALLOCATE(WALL_DUMMY)
      M%N_WALL_CELLS_DIM = M%N_WALL_CELLS_DIM + N_NEW_STORAGE_SLOTS
   ENDIF

ELSEIF (PRESENT(CFACE_INDEX)) THEN

   OS => MESHES(NM)%CFACE_STORAGE(SURF_INDEX)

   ! Check if there are enough CFACEs

   IF (CFACE_INDEX>M%N_CFACE_CELLS_DIM) THEN
      ALLOCATE(CFACE_DUMMY(0:M%N_CFACE_CELLS_DIM))
      CFACE_DUMMY = M%CFACE
      DEALLOCATE(M%CFACE)
      ALLOCATE(M%CFACE(0:M%N_CFACE_CELLS_DIM+N_NEW_STORAGE_SLOTS))
      M%CFACE(0:MESHES(NM)%N_CFACE_CELLS_DIM) = CFACE_DUMMY
      DEALLOCATE(CFACE_DUMMY)
      M%N_CFACE_CELLS_DIM = M%N_CFACE_CELLS_DIM + N_NEW_STORAGE_SLOTS
   ENDIF

ELSEIF (PRESENT(LP_INDEX)) THEN

   LPC => LAGRANGIAN_PARTICLE_CLASS(LPC_INDEX)
   OS  => MESHES(NM)%PARTICLE_STORAGE(LPC_INDEX)

   ! Check if there are enough particles

   IF (LP_INDEX>MESHES(NM)%NLPDIM) THEN
      ALLOCATE(LP_DUMMY(MESHES(NM)%NLPDIM))
      LP_DUMMY = MESHES(NM)%LAGRANGIAN_PARTICLE
      DEALLOCATE(MESHES(NM)%LAGRANGIAN_PARTICLE)
      ALLOCATE(MESHES(NM)%LAGRANGIAN_PARTICLE(MESHES(NM)%NLPDIM+N_NEW_STORAGE_SLOTS))
      MESHES(NM)%LAGRANGIAN_PARTICLE(1:MESHES(NM)%NLPDIM) = LP_DUMMY
      DEALLOCATE(LP_DUMMY)
      MESHES(NM)%NLPDIM = MESHES(NM)%NLPDIM + N_NEW_STORAGE_SLOTS
   ENDIF

ELSE

   RETURN

ENDIF

! Determine if this is to be a new entry or old

STORAGE_INDEX_OLD  = 0
STORAGE_INDEX_OPEN = HUGE(STORAGE_INDEX_OPEN)

! If there is a specified particle TAG, find its storage slot

IF (PRESENT(TAG)) THEN
   IF (PRESENT(NEW_TAG) .AND. .NOT.NEW_TAG) THEN
      DO I=1,OS%N_STORAGE_SLOTS
         IF (OS%INTEGERS(1,I)==TAG) THEN
            STORAGE_INDEX_OLD = I
            EXIT
         ENDIF
      ENDDO
   ENDIF
ENDIF

! If there is an existing storage slot for the wall cell, use it

IF (PRESENT(WALL_INDEX) .AND. .NOT.INITIALIZATION_PHASE) THEN
   DO I=1,OS%N_STORAGE_SLOTS
      IF (OS%INTEGERS(1,I)==WALL_INDEX) THEN
         STORAGE_INDEX_OLD = I
         EXIT
      ENDIF
   ENDDO
ENDIF

! Look for the first available storage slot if it fails to find an existing one

IF (STORAGE_INDEX_OLD==0) THEN
   NAS = OS%NEXT_AVAILABLE_SLOT
   DO I=NAS,OS%N_STORAGE_SLOTS
      IF (OS%INTEGERS(1,I)==0) THEN
         STORAGE_INDEX_OPEN = MIN(I,STORAGE_INDEX_OPEN)
         OS%NEXT_AVAILABLE_SLOT = I+1
         EXIT
      ENDIF
   ENDDO
ENDIF

! Decide which storage index to use

IF (STORAGE_INDEX_OLD>0) THEN
   STORAGE_INDEX = STORAGE_INDEX_OLD
   NEW_SLOT = .FALSE.
ELSEIF (STORAGE_INDEX_OPEN<HUGE(STORAGE_INDEX_OPEN)) THEN
   STORAGE_INDEX = STORAGE_INDEX_OPEN
   NEW_SLOT = .TRUE.
ELSE
   STORAGE_INDEX = 0
   NEW_SLOT = .TRUE.
ENDIF

! Re-allocate storage arrays if necessary

IF (STORAGE_INDEX==0) THEN
   N_STORAGE_SLOTS_OLD = OS%N_STORAGE_SLOTS
   IF (PRESENT(LP_INDEX))    CALL REALLOCATE_STORAGE_ARRAYS(NM,1,LPC_INDEX,N_NEW_STORAGE_SLOTS)
   IF (PRESENT(WALL_INDEX))  CALL REALLOCATE_STORAGE_ARRAYS(NM,4,SURF_INDEX,N_NEW_STORAGE_SLOTS)
   IF (PRESENT(CFACE_INDEX)) CALL REALLOCATE_STORAGE_ARRAYS(NM,5,SURF_INDEX,N_NEW_STORAGE_SLOTS)
   DO I=1,N_STORAGE_SLOTS_OLD
      IF (PRESENT(WALL_INDEX))  OLD_INDEX = OS%INTEGERS(1,I)
      IF (PRESENT(CFACE_INDEX)) OLD_INDEX = OS%INTEGERS(1,I)
      IF (PRESENT(LP_INDEX))    OLD_INDEX = OS%INTEGERS(2,I)
      STORAGE_INDEX = I
      IF (PRESENT(WALL_INDEX)  .AND. OLD_INDEX>0) CALL WALL_POINTERS(OLD_INDEX,.FALSE.)
      IF (PRESENT(CFACE_INDEX) .AND. OLD_INDEX>0) CALL CFACE_POINTERS(OLD_INDEX,.FALSE.)
      IF (PRESENT(LP_INDEX)    .AND. OLD_INDEX>0) CALL PARTICLE_POINTERS(OLD_INDEX,.FALSE.)
   ENDDO
   STORAGE_INDEX = N_STORAGE_SLOTS_OLD + 1
ENDIF

! Assign values to the pointer arrays

IF (PRESENT(WALL_INDEX))  CALL WALL_POINTERS(WALL_INDEX,NEW_SLOT)
IF (PRESENT(CFACE_INDEX)) CALL CFACE_POINTERS(CFACE_INDEX,NEW_SLOT)
IF (PRESENT(LP_INDEX))    CALL PARTICLE_POINTERS(LP_INDEX,NEW_SLOT)

CONTAINS


SUBROUTINE PARTICLE_POINTERS(LP_INDEX_LOCAL,NEW)

INTEGER, INTENT(IN) :: LP_INDEX_LOCAL
LOGICAL, INTENT(IN) :: NEW

LP => MESHES(NM)%LAGRANGIAN_PARTICLE(LP_INDEX_LOCAL)

! Assign integer pointers and values

LP%TAG                 => OS%INTEGERS( 1,STORAGE_INDEX) ; IF (NEW) LP%TAG = TAG
LP%ARRAY_INDEX         => OS%INTEGERS( 2,STORAGE_INDEX) ;          LP%ARRAY_INDEX = LP_INDEX_LOCAL
LP%STORAGE_INDEX       => OS%INTEGERS( 3,STORAGE_INDEX) ; IF (NEW) LP%STORAGE_INDEX = STORAGE_INDEX
LP%CLASS_INDEX         => OS%INTEGERS( 4,STORAGE_INDEX) ; IF (NEW) LP%CLASS_INDEX = LPC_INDEX
LP%ORIENTATION_INDEX   => OS%INTEGERS( 5,STORAGE_INDEX) ; IF (NEW) LP%ORIENTATION_INDEX = 0
LP%WALL_INDEX          => OS%INTEGERS( 6,STORAGE_INDEX) ; IF (NEW) LP%WALL_INDEX = 0
LP%FACE_INDEX          => OS%INTEGERS( 7,STORAGE_INDEX) ; IF (NEW) LP%FACE_INDEX = 0
LP%DUCT_INDEX          => OS%INTEGERS( 8,STORAGE_INDEX) ; IF (NEW) LP%DUCT_INDEX = 0
LP%DUCT_CELL_INDEX     => OS%INTEGERS( 9,STORAGE_INDEX) ; IF (NEW) LP%DUCT_CELL_INDEX = 0

! Assign logical pointers and values

LP%SHOW              => OS%LOGICALS( 1,STORAGE_INDEX) ; IF (NEW) LP%SHOW  = .FALSE.
LP%SPLAT             => OS%LOGICALS( 2,STORAGE_INDEX) ; IF (NEW) LP%SPLAT = .FALSE.
LP%EMBER             => OS%LOGICALS( 3,STORAGE_INDEX) ; IF (NEW) LP%EMBER = .FALSE.

! Assign real pointers and values

LP%X                 => OS%REALS( 1,STORAGE_INDEX) ; IF (NEW) LP%X       = 0._EB
LP%Y                 => OS%REALS( 2,STORAGE_INDEX) ; IF (NEW) LP%Y       = 0._EB
LP%Z                 => OS%REALS( 3,STORAGE_INDEX) ; IF (NEW) LP%Z       = 0._EB
LP%U                 => OS%REALS( 4,STORAGE_INDEX) ; IF (NEW) LP%U       = 0._EB
LP%V                 => OS%REALS( 5,STORAGE_INDEX) ; IF (NEW) LP%V       = 0._EB
LP%W                 => OS%REALS( 6,STORAGE_INDEX) ; IF (NEW) LP%W       = 0._EB
LP%PWT               => OS%REALS( 7,STORAGE_INDEX) ; IF (NEW) LP%PWT     = 1._EB
LP%ACCEL_X           => OS%REALS( 8,STORAGE_INDEX) ; IF (NEW) LP%ACCEL_X = 0._EB
LP%ACCEL_Y           => OS%REALS( 9,STORAGE_INDEX) ; IF (NEW) LP%ACCEL_Y = 0._EB
LP%ACCEL_Z           => OS%REALS(10,STORAGE_INDEX) ; IF (NEW) LP%ACCEL_Z = 0._EB
LP%RE                => OS%REALS(11,STORAGE_INDEX) ; IF (NEW) LP%RE      = 0._EB
LP%MASS              => OS%REALS(12,STORAGE_INDEX) ; IF (NEW) LP%MASS    = 0._EB
LP%T_INSERT          => OS%REALS(13,STORAGE_INDEX) ; IF (NEW) LP%T_INSERT= 0._EB
LP%DX                => OS%REALS(14,STORAGE_INDEX) ; IF (NEW) LP%DX      = 0._EB
LP%DY                => OS%REALS(15,STORAGE_INDEX) ; IF (NEW) LP%DY      = 0._EB
LP%DZ                => OS%REALS(16,STORAGE_INDEX) ; IF (NEW) LP%DZ      = 0._EB
LP%M_DOT             => OS%REALS(17,STORAGE_INDEX) ; IF (NEW) LP%M_DOT   = 0._EB

CALL ONE_D_POINTERS(LP_INDEX_LOCAL,NEW,N_PARTICLE_SCALAR_REALS,N_PARTICLE_SCALAR_INTEGERS,N_PARTICLE_SCALAR_LOGICALS)

END SUBROUTINE PARTICLE_POINTERS


SUBROUTINE WALL_POINTERS(WALL_INDEX_LOCAL,NEW)

INTEGER, INTENT(IN) :: WALL_INDEX_LOCAL
LOGICAL, INTENT(IN) :: NEW
TYPE (WALL_TYPE), POINTER :: WC

WC => MESHES(NM)%WALL(WALL_INDEX_LOCAL)

! Assign and initialize integers

WC%WALL_INDEX         => OS%INTEGERS( 1,STORAGE_INDEX) ; IF (NEW) WC%WALL_INDEX = WALL_INDEX
WC%SURF_INDEX         => OS%INTEGERS( 2,STORAGE_INDEX) ; IF (NEW) WC%SURF_INDEX = SURF_INDEX
WC%BACK_INDEX         => OS%INTEGERS( 3,STORAGE_INDEX) ; IF (NEW) WC%BACK_INDEX = 0
WC%BACK_MESH          => OS%INTEGERS( 4,STORAGE_INDEX) ; IF (NEW) WC%BACK_MESH  = NM
WC%BOUNDARY_TYPE      => OS%INTEGERS( 5,STORAGE_INDEX) ; IF (NEW) WC%BOUNDARY_TYPE = NULL_BOUNDARY
WC%SURF_INDEX_ORIG    => OS%INTEGERS( 6,STORAGE_INDEX) ; IF (NEW) WC%SURF_INDEX_ORIG = 0
WC%OBST_INDEX         => OS%INTEGERS( 7,STORAGE_INDEX) ; IF (NEW) WC%OBST_INDEX = 0
WC%PRESSURE_BC_INDEX  => OS%INTEGERS( 8,STORAGE_INDEX) ; IF (NEW) WC%PRESSURE_BC_INDEX = DIRICHLET
WC%LAPLACE_BC_INDEX   => OS%INTEGERS( 9,STORAGE_INDEX) ; IF (NEW) WC%LAPLACE_BC_INDEX = NEUMANN
WC%VENT_INDEX         => OS%INTEGERS(10,STORAGE_INDEX) ; IF (NEW) WC%VENT_INDEX = 0
WC%NODE_INDEX         => OS%INTEGERS(11,STORAGE_INDEX) ; IF (NEW) WC%NODE_INDEX = 0
WC%JD11_INDEX         => OS%INTEGERS(12,STORAGE_INDEX) ; IF (NEW) WC%JD11_INDEX = 0
WC%JD12_INDEX         => OS%INTEGERS(13,STORAGE_INDEX) ; IF (NEW) WC%JD12_INDEX = 0
WC%JD21_INDEX         => OS%INTEGERS(14,STORAGE_INDEX) ; IF (NEW) WC%JD21_INDEX = 0
WC%JD22_INDEX         => OS%INTEGERS(15,STORAGE_INDEX) ; IF (NEW) WC%JD22_INDEX = 0

! Assign and initialize reals

WC%DUWDT              => OS%REALS( 1,STORAGE_INDEX) ; IF (NEW) WC%DUWDT = 0._EB
WC%EW                 => OS%REALS( 2,STORAGE_INDEX) ; IF (NEW) WC%EW = 0._EB
WC%UW0                => OS%REALS( 3,STORAGE_INDEX) ; IF (NEW) WC%UW0 = 0._EB
WC%X                  => OS%REALS( 4,STORAGE_INDEX)
WC%Y                  => OS%REALS( 5,STORAGE_INDEX)
WC%Z                  => OS%REALS( 6,STORAGE_INDEX)
WC%VEG_HEIGHT         => OS%REALS( 7,STORAGE_INDEX)
WC%VEL_ERR_NEW        => OS%REALS( 8,STORAGE_INDEX) ; IF (NEW) WC%VEL_ERR_NEW = 0._EB
WC%V_DEP              => OS%REALS( 9,STORAGE_INDEX) ; IF (NEW) WC%V_DEP = 0._EB
WC%Q_LEAK             => OS%REALS(10,STORAGE_INDEX) ; IF (NEW) WC%Q_LEAK = 0._EB

I1 = N_WALL_SCALAR_REALS+1 ; I2 = I1 + SF%WALL_REALS_ARRAY_SIZE(1) - 1  ; WC%A_LP_MPUA(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
I1 = I2+1 ; I2 = I1 + SF%WALL_REALS_ARRAY_SIZE(2) - 1  ; WC%LP_CPUA(1:I2-I1+1)         => OS%REALS(I1:I2,STORAGE_INDEX)
I1 = I2+1 ; I2 = I1 + SF%WALL_REALS_ARRAY_SIZE(3) - 1  ; WC%LP_MPUA(1:I2-I1+1)         => OS%REALS(I1:I2,STORAGE_INDEX)
I1 = I2+1 ; I2 = I1 + SF%WALL_REALS_ARRAY_SIZE(4) - 1  ; WC%AWM_AEROSOL(1:I2-I1+1)     => OS%REALS(I1:I2,STORAGE_INDEX)
I1 = I2+1 ; I2 = I1 + SF%WALL_REALS_ARRAY_SIZE(5) - 1  ; WC%VEG_FUELMASS_L(1:I2-I1+1)  => OS%REALS(I1:I2,STORAGE_INDEX)
I1 = I2+1 ; I2 = I1 + SF%WALL_REALS_ARRAY_SIZE(6) - 1  ; WC%VEG_MOISTMASS_L(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
I1 = I2+1 ; I2 = I1 + SF%WALL_REALS_ARRAY_SIZE(7) - 1  ; WC%VEG_TMP_L(0:I2-I1)         => OS%REALS(I1:I2,STORAGE_INDEX)

CALL ONE_D_POINTERS(WALL_INDEX_LOCAL,NEW,N_WALL_SCALAR_REALS+SUM(SF%WALL_REALS_ARRAY_SIZE(1:7)), &
                    N_WALL_SCALAR_INTEGERS,N_WALL_SCALAR_LOGICALS)

END SUBROUTINE WALL_POINTERS


SUBROUTINE CFACE_POINTERS(CFACE_INDEX_LOCAL,NEW)

INTEGER, INTENT(IN) :: CFACE_INDEX_LOCAL
LOGICAL, INTENT(IN) :: NEW
TYPE (CFACE_TYPE), POINTER :: CFA

CFA => MESHES(NM)%CFACE(CFACE_INDEX_LOCAL)

! Assign and initialize integers

CFA%CFACE_INDEX        => OS%INTEGERS(1,STORAGE_INDEX) ; IF (NEW) CFA%CFACE_INDEX   = CFACE_INDEX
CFA%SURF_INDEX         => OS%INTEGERS(2,STORAGE_INDEX) ; IF (NEW) CFA%SURF_INDEX    = SURF_INDEX
CFA%BOUNDARY_TYPE      => OS%INTEGERS(3,STORAGE_INDEX) ; IF (NEW) CFA%BOUNDARY_TYPE = NULL_BOUNDARY
CFA%CUT_FACE_IND1      => OS%INTEGERS(4,STORAGE_INDEX) ; IF (NEW) CFA%CUT_FACE_IND1 = -11
CFA%CUT_FACE_IND2      => OS%INTEGERS(5,STORAGE_INDEX) ; IF (NEW) CFA%CUT_FACE_IND2 = -11

! Assign and initialize reals

CFA%AREA               => OS%REALS(1,STORAGE_INDEX)
CFA%X                  => OS%REALS(2,STORAGE_INDEX)
CFA%Y                  => OS%REALS(3,STORAGE_INDEX)
CFA%Z                  => OS%REALS(4,STORAGE_INDEX)
CFA%NVEC(1:3)          => OS%REALS(5:7,STORAGE_INDEX)
CFA%VEL_ERR_NEW        => OS%REALS(8,STORAGE_INDEX)  ; IF (NEW) CFA%VEL_ERR_NEW = 0._EB
CFA%UW0                => OS%REALS(9,STORAGE_INDEX)  ; IF (NEW) CFA%UW0         = 0._EB
CFA%V_DEP              => OS%REALS(10,STORAGE_INDEX) ; IF (NEW) CFA%V_DEP       = 0._EB

I1=N_CFACE_SCALAR_REALS+1; I2=I1+SF%CFACE_REALS_ARRAY_SIZE(1)-1;
!                                                CFA%A_LP_MPUA(1:I2-I1+1)       => OS%REALS(I1:I2,STORAGE_INDEX)
!I1=I2+1 ; I2=I1+SF%CFACE_REALS_ARRAY_SIZE(2)-1; CFA%LP_CPUA(1:I2-I1+1)         => OS%REALS(I1:I2,STORAGE_INDEX)
!I1=I2+1 ; I2=I1+SF%CFACE_REALS_ARRAY_SIZE(3)-1; CFA%LP_MPUA(1:I2-I1+1)         => OS%REALS(I1:I2,STORAGE_INDEX)
I1=I2+1 ; I2=I1+SF%CFACE_REALS_ARRAY_SIZE(4)-1; CFA%AWM_AEROSOL(1:I2-I1+1)     => OS%REALS(I1:I2,STORAGE_INDEX)
!I1=I2+1 ; I2=I1+SF%CFACE_REALS_ARRAY_SIZE(5)-1; CFA%VEG_FUELMASS_L(1:I2-I1+1)  => OS%REALS(I1:I2,STORAGE_INDEX)
!I1=I2+1 ; I2=I1+SF%CFACE_REALS_ARRAY_SIZE(6)-1; CFA%VEG_MOISTMASS_L(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
!I1=I2+1 ; I2=I1+SF%CFACE_REALS_ARRAY_SIZE(7)-1; CFA%VEG_TMP_L(0:I2-I1)         => OS%REALS(I1:I2,STORAGE_INDEX)


CALL ONE_D_POINTERS(CFACE_INDEX_LOCAL,NEW,N_CFACE_SCALAR_REALS+SUM(SF%CFACE_REALS_ARRAY_SIZE(1:7)),&
                    N_CFACE_SCALAR_INTEGERS,N_CFACE_SCALAR_LOGICALS)

END SUBROUTINE CFACE_POINTERS


SUBROUTINE ONE_D_POINTERS(ONE_D_INDEX_LOCAL,NEW,RC,IC,LC)

USE GLOBAL_CONSTANTS, ONLY: RADIATION,TMPM,RHOA,RSUM0
INTEGER, INTENT(IN) :: ONE_D_INDEX_LOCAL,RC,IC,LC
LOGICAL, INTENT(IN) :: NEW
INTEGER :: N,NN
REAL(EB) :: X_POS
TYPE(ONE_D_M_AND_E_XFER_TYPE), POINTER :: ONE_D

IF (PRESENT(WALL_INDEX)) THEN
   ONE_D=>MESHES(NM)%WALL(ONE_D_INDEX_LOCAL)%ONE_D
ELSEIF (PRESENT(CFACE_INDEX)) THEN
   ONE_D=>MESHES(NM)%CFACE(ONE_D_INDEX_LOCAL)%ONE_D
ELSEIF(PRESENT(LP_INDEX)) THEN
   ONE_D=>MESHES(NM)%LAGRANGIAN_PARTICLE(ONE_D_INDEX_LOCAL)%ONE_D
ENDIF

! Assign and initialize integers

ONE_D%ARRAY_INDEX     => OS%INTEGERS(IC+ 1,STORAGE_INDEX) ; IF (NEW) ONE_D%ARRAY_INDEX = ONE_D_INDEX_LOCAL
ONE_D%STORAGE_INDEX   => OS%INTEGERS(IC+ 2,STORAGE_INDEX) ; IF (NEW) ONE_D%STORAGE_INDEX = STORAGE_INDEX
ONE_D%II              => OS%INTEGERS(IC+ 3,STORAGE_INDEX) ; IF (NEW) ONE_D%II            = 0
ONE_D%JJ              => OS%INTEGERS(IC+ 4,STORAGE_INDEX) ; IF (NEW) ONE_D%JJ            = 0
ONE_D%KK              => OS%INTEGERS(IC+ 5,STORAGE_INDEX) ; IF (NEW) ONE_D%KK            = 0
ONE_D%IIG             => OS%INTEGERS(IC+ 6,STORAGE_INDEX) ; IF (NEW) ONE_D%IIG           = 0
ONE_D%JJG             => OS%INTEGERS(IC+ 7,STORAGE_INDEX) ; IF (NEW) ONE_D%JJG           = 0
ONE_D%KKG             => OS%INTEGERS(IC+ 8,STORAGE_INDEX) ; IF (NEW) ONE_D%KKG           = 0
ONE_D%IOR             => OS%INTEGERS(IC+ 9,STORAGE_INDEX) ; IF (NEW) ONE_D%IOR           = 0
ONE_D%PRESSURE_ZONE   => OS%INTEGERS(IC+10,STORAGE_INDEX) ; IF (NEW) ONE_D%PRESSURE_ZONE = 0

I1 = IC+1+N_ONE_D_SCALAR_INTEGERS ; I2 = I1 + SF%ONE_D_INTEGERS_ARRAY_SIZE(1) - 1
ONE_D%N_LAYER_CELLS(1:I2-I1+1) => OS%INTEGERS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%N_LAYER_CELLS(1:SF%N_LAYERS) = SF%N_LAYER_CELLS(1:SF%N_LAYERS)

! Assign and initialize logicals

ONE_D%BURNAWAY        => OS%LOGICALS(LC+1,STORAGE_INDEX) ; IF (NEW) ONE_D%BURNAWAY = .FALSE.
ONE_D%CHANGE_THICKNESS => OS%LOGICALS(LC+2,STORAGE_INDEX) ; IF (NEW) ONE_D%CHANGE_THICKNESS = .FALSE.

! Assign and initialize reals

ONE_D%AREA            => OS%REALS(RC+ 1,STORAGE_INDEX) ; IF (NEW) ONE_D%AREA            = 0._EB
ONE_D%HEAT_TRANS_COEF => OS%REALS(RC+ 2,STORAGE_INDEX) ; IF (NEW) ONE_D%HEAT_TRANS_COEF = 0._EB
ONE_D%QCONF           => OS%REALS(RC+ 3,STORAGE_INDEX) ; IF (NEW) ONE_D%QCONF           = 0._EB
IF (RADIATION) THEN
   ONE_D%QRADIN          => OS%REALS(RC+ 4,STORAGE_INDEX) ; IF (NEW) ONE_D%QRADIN          = SF%EMISSIVITY*SIGMA*TMPA4
   ONE_D%QRADOUT         => OS%REALS(RC+ 5,STORAGE_INDEX) ; IF (NEW) ONE_D%QRADOUT         = SF%EMISSIVITY*SIGMA*TMPA4
ELSE
   ONE_D%QRADIN          => OS%REALS(RC+ 4,STORAGE_INDEX) ; IF (NEW) ONE_D%QRADIN          = 0._EB
   ONE_D%QRADOUT         => OS%REALS(RC+ 5,STORAGE_INDEX) ; IF (NEW) ONE_D%QRADOUT         = 0._EB
ENDIF
ONE_D%EMISSIVITY      => OS%REALS(RC+ 6,STORAGE_INDEX) ; IF (NEW) ONE_D%EMISSIVITY      = SF%EMISSIVITY
ONE_D%AREA_ADJUST     => OS%REALS(RC+ 7,STORAGE_INDEX) ; IF (NEW) ONE_D%AREA_ADJUST     = 1._EB
ONE_D%T_IGN           => OS%REALS(RC+ 8,STORAGE_INDEX) ; IF (NEW) ONE_D%T_IGN           = T_BEGIN
ONE_D%TMP_F           => OS%REALS(RC+ 9,STORAGE_INDEX) ; IF (NEW) ONE_D%TMP_F           = SF%TMP_FRONT
ONE_D%TMP_B           => OS%REALS(RC+10,STORAGE_INDEX) ; IF (NEW) ONE_D%TMP_B           = SF%TMP_INNER(SF%N_LAYERS)
ONE_D%UW              => OS%REALS(RC+11,STORAGE_INDEX) ; IF (NEW) ONE_D%UW              = 0._EB
ONE_D%UWS             => OS%REALS(RC+12,STORAGE_INDEX) ; IF (NEW) ONE_D%UWS             = 0._EB
ONE_D%RSUM_G          => OS%REALS(RC+13,STORAGE_INDEX) ; IF (NEW) ONE_D%RSUM_G          = RSUM0
ONE_D%TMP_G           => OS%REALS(RC+14,STORAGE_INDEX) ; IF (NEW) ONE_D%TMP_G           = TMPA
ONE_D%RHO_G           => OS%REALS(RC+15,STORAGE_INDEX) ; IF (NEW) ONE_D%RHO_G           = RHOA
ONE_D%U_TANG          => OS%REALS(RC+16,STORAGE_INDEX) ; IF (NEW) ONE_D%U_TANG          = 0._EB
ONE_D%RHO_F           => OS%REALS(RC+17,STORAGE_INDEX) ; IF (NEW) ONE_D%RHO_F           = RHOA
ONE_D%RDN             => OS%REALS(RC+18,STORAGE_INDEX) ; IF (NEW) ONE_D%RDN             = 1._EB
ONE_D%MU_G            => OS%REALS(RC+19,STORAGE_INDEX) ; IF (NEW) ONE_D%MU_G            = 0.1_EB
ONE_D%K_G             => OS%REALS(RC+20,STORAGE_INDEX) ; IF (NEW) ONE_D%K_G             = 0.1_EB
ONE_D%U_TAU           => OS%REALS(RC+21,STORAGE_INDEX) ; IF (NEW) ONE_D%U_TAU           = 0._EB
ONE_D%Y_PLUS          => OS%REALS(RC+22,STORAGE_INDEX) ; IF (NEW) ONE_D%Y_PLUS          = 1._EB
ONE_D%Z_STAR          => OS%REALS(RC+23,STORAGE_INDEX) ; IF (NEW) ONE_D%Z_STAR          = 1._EB

I1 = RC+1+N_ONE_D_SCALAR_REALS ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(1) - 1
ONE_D%MASSFLUX_SPEC(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%MASSFLUX_SPEC(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(2) - 1
ONE_D%MASSFLUX_MATL(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%MASSFLUX_MATL(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(3) - 1
ONE_D%MASSFLUX(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%MASSFLUX(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(4) - 1
ONE_D%X(0:I2-I1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%X(0:SF%N_CELLS_INI) = SF%X_S(0:SF%N_CELLS_INI)

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(5) - 1
ONE_D%TMP(0:I2-I1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%TMP(0:SF%N_CELLS_INI+1) = SF%TMP_INNER(1)

IF (NEW) THEN
   IF (SF%THERMALLY_THICK) THEN
      IF (SF%RAMP_T_I_INDEX > 0) THEN
         DO N=0,SF%N_CELLS_INI
            X_POS = MAX(0._EB,MIN(RAMPS(SF%RAMP_T_I_INDEX)%SPAN,SF%X_S(N)-RAMPS(SF%RAMP_T_I_INDEX)%T_MIN))
            ONE_D%TMP(N)=RAMPS(SF%RAMP_T_I_INDEX)%INTERPOLATED_DATA(NINT(X_POS*RAMPS(SF%RAMP_T_I_INDEX)%RDT)) + TMPM
         ENDDO
         ONE_D%TMP(SF%N_CELLS_INI+1) = ONE_D%TMP(SF%N_CELLS_INI)
         ONE_D%TMP_F = ONE_D%TMP(0)
         ONE_D%TMP_B = ONE_D%TMP(SF%N_CELLS_INI+1)
      ELSE
         ONE_D%TMP(0) = SF%TMP_INNER(1)
         ONE_D%TMP(SF%N_CELLS_INI+1) = SF%TMP_INNER(SF%N_LAYERS)
         NN=1
         DO N=1,SF%N_LAYERS
            ONE_D%TMP(NN:NN+SF%N_LAYER_CELLS(N)-1) = SF%TMP_INNER(N)
            NN = NN + SF%N_LAYER_CELLS(N)
         ENDDO
      ENDIF
   ELSE
      ONE_D%TMP(0:SF%N_CELLS_INI+1) = SF%TMP_INNER(1)
   ENDIF
ENDIF

IF (ALLOCATED(ONE_D%MATL_COMP)) DEALLOCATE(ONE_D%MATL_COMP)
ALLOCATE(ONE_D%MATL_COMP(1:SF%N_MATL))
DO NN=1,SF%N_MATL
   I1 = I2+1 ; I2 = I1 + (SF%N_CELLS_MAX+2) - 1
   ONE_D%MATL_COMP(NN)%RHO(0:SF%N_CELLS_MAX+1) => OS%REALS(I1:I2,STORAGE_INDEX)
   IF (NEW) ONE_D%MATL_COMP(NN)%RHO(0:SF%N_CELLS_INI+1) = SF%RHO_0(0:SF%N_CELLS_INI+1,NN)
ENDDO

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(7) - 1
ONE_D%LAYER_THICKNESS(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%LAYER_THICKNESS(1:SF%N_LAYERS) = SF%LAYER_THICKNESS(1:SF%N_LAYERS)

IF (ALLOCATED(ONE_D%BAND)) DEALLOCATE(ONE_D%BAND)
ALLOCATE(ONE_D%BAND(1:SF%NSB))
DO NN=1,SF%NSB
   I1 = I2+1 ; I2 = I1 + SF%NRA - 1
   ONE_D%BAND(NN)%ILW(1:SF%NRA) => OS%REALS(I1:I2,STORAGE_INDEX)
   IF (NEW) ONE_D%BAND(NN)%ILW(1:SF%NRA) = 0._EB
ENDDO

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(9) - 1
ONE_D%IL(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%IL(1:SF%NSB) = SIGMA*TMPA4/PI

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(10) - 1
ONE_D%ZZ_G(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%ZZ_G(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(11) - 1
ONE_D%ZZ_F(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%ZZ_F(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(12) - 1
ONE_D%RHO_D_F(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%RHO_D_F(1:I2-I1+1) = 0._EB

I1 = I2+1 ; I2 = I1 + SF%ONE_D_REALS_ARRAY_SIZE(13) - 1
ONE_D%RHO_D_DZDN_F(1:I2-I1+1) => OS%REALS(I1:I2,STORAGE_INDEX)
IF (NEW) ONE_D%RHO_D_DZDN_F(1:I2-I1+1) = 0._EB

END SUBROUTINE ONE_D_POINTERS

END SUBROUTINE ALLOCATE_STORAGE


SUBROUTINE REALLOCATE_STORAGE_ARRAYS(NM,MODE,INDX,N_NEW_STORAGE_SLOTS,NOM)

INTEGER, INTENT(IN) :: MODE,NM,INDX,N_NEW_STORAGE_SLOTS
INTEGER, INTENT(IN), OPTIONAL :: NOM
INTEGER :: REALS_DIMENSION,INTEGERS_DIMENSION,LOGICALS_DIMENSION
TYPE (STORAGE_TYPE), POINTER :: OS
TYPE (SURFACE_TYPE), POINTER :: SF
TYPE (LAGRANGIAN_PARTICLE_CLASS_TYPE), POINTER :: LPC
REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: DUMMY_REALS
INTEGER, ALLOCATABLE, DIMENSION(:,:) :: DUMMY_INTEGERS
LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: DUMMY_LOGICALS

SELECT CASE(MODE)
   CASE(1:3)
      LPC => LAGRANGIAN_PARTICLE_CLASS(INDX)
      IF (MODE==1) OS  => MESHES(NM)%PARTICLE_STORAGE(INDX)
      IF (MODE==2) OS  => MESHES(NM)%OMESH(NOM)%ORPHAN_PARTICLE_STORAGE(INDX)
      IF (MODE==3) OS  => MESHES(NM)%OMESH(NOM)%ADOPT_PARTICLE_STORAGE(INDX)
      REALS_DIMENSION    = LPC%N_STORAGE_REALS
      INTEGERS_DIMENSION = LPC%N_STORAGE_INTEGERS
      LOGICALS_DIMENSION = LPC%N_STORAGE_LOGICALS
   CASE(4)
      SF => SURFACE(INDX)
      OS => MESHES(NM)%WALL_STORAGE(INDX)
      REALS_DIMENSION    = SF%N_WALL_STORAGE_REALS
      INTEGERS_DIMENSION = SF%N_WALL_STORAGE_INTEGERS
      LOGICALS_DIMENSION = SF%N_WALL_STORAGE_LOGICALS
   CASE(5)
      SF => SURFACE(INDX)
      OS => MESHES(NM)%CFACE_STORAGE(INDX)
      REALS_DIMENSION    = SF%N_CFACE_STORAGE_REALS
      INTEGERS_DIMENSION = SF%N_CFACE_STORAGE_INTEGERS
      LOGICALS_DIMENSION = SF%N_CFACE_STORAGE_LOGICALS
END SELECT

IF (OS%N_STORAGE_SLOTS>0) THEN
   ALLOCATE(DUMMY_REALS(REALS_DIMENSION,OS%N_STORAGE_SLOTS))
   ALLOCATE(DUMMY_INTEGERS(INTEGERS_DIMENSION,OS%N_STORAGE_SLOTS))
   ALLOCATE(DUMMY_LOGICALS(LOGICALS_DIMENSION,OS%N_STORAGE_SLOTS))
   DUMMY_REALS    = OS%REALS
   DUMMY_INTEGERS = OS%INTEGERS
   DUMMY_LOGICALS = OS%LOGICALS
   DEALLOCATE(OS%REALS)
   DEALLOCATE(OS%INTEGERS)
   DEALLOCATE(OS%LOGICALS)
ENDIF

ALLOCATE(OS%REALS(REALS_DIMENSION,OS%N_STORAGE_SLOTS+N_NEW_STORAGE_SLOTS))
ALLOCATE(OS%INTEGERS(INTEGERS_DIMENSION,OS%N_STORAGE_SLOTS+N_NEW_STORAGE_SLOTS))
ALLOCATE(OS%LOGICALS(LOGICALS_DIMENSION,OS%N_STORAGE_SLOTS+N_NEW_STORAGE_SLOTS))

OS%REALS    = 0._EB
OS%INTEGERS = 0
OS%LOGICALS = .FALSE.

IF (OS%N_STORAGE_SLOTS>0) THEN
   OS%REALS(1:REALS_DIMENSION,1:OS%N_STORAGE_SLOTS) = DUMMY_REALS
   OS%INTEGERS(1:INTEGERS_DIMENSION,1:OS%N_STORAGE_SLOTS) = DUMMY_INTEGERS
   OS%LOGICALS(1:LOGICALS_DIMENSION,1:OS%N_STORAGE_SLOTS) = DUMMY_LOGICALS
ENDIF

IF (OS%N_STORAGE_SLOTS>0) THEN
   DEALLOCATE(DUMMY_REALS)
   DEALLOCATE(DUMMY_INTEGERS)
   DEALLOCATE(DUMMY_LOGICALS)
ENDIF

OS%N_STORAGE_SLOTS = OS%N_STORAGE_SLOTS + N_NEW_STORAGE_SLOTS

END SUBROUTINE REALLOCATE_STORAGE_ARRAYS


SUBROUTINE GET_LAGRANGIAN_PARTICLE_INDEX(NM,LPC_INDEX,LP_TAG,LP_INDEX)

! On mesh NM, find the index of the particle with the given TAG and LPC (Class). If not found, return 0.

USE GLOBAL_CONSTANTS, ONLY: EVACUATION_ONLY
INTEGER, INTENT(IN) :: NM,LPC_INDEX,LP_TAG
INTEGER, INTENT(OUT) :: LP_INDEX
INTEGER :: I

LP_INDEX = 0
IF (EVACUATION_ONLY(NM)) RETURN
DO I=1,MESHES(NM)%PARTICLE_STORAGE(LPC_INDEX)%N_STORAGE_SLOTS
   IF (MESHES(NM)%PARTICLE_STORAGE(LPC_INDEX)%INTEGERS(1,I)==LP_TAG) THEN
      LP_INDEX = MESHES(NM)%PARTICLE_STORAGE(LPC_INDEX)%INTEGERS(2,I)
      EXIT
   ENDIF
ENDDO

END SUBROUTINE GET_LAGRANGIAN_PARTICLE_INDEX


END MODULE MEMORY_FUNCTIONS



MODULE GEOMETRY_FUNCTIONS

! Functions for manipulating geometry

USE PRECISION_PARAMETERS
USE MESH_VARIABLES
USE GLOBAL_CONSTANTS
IMPLICIT NONE

CONTAINS


SUBROUTINE SEARCH_OTHER_MESHES(XX,YY,ZZ,NOM,IIO,JJO,KKO)

! Given the point (XX,YY,ZZ), determine which other mesh it intersects and what its indices are.

REAL(EB), INTENT(IN) :: XX,YY,ZZ
REAL(EB) :: XI,YJ,ZK
INTEGER, INTENT(OUT) :: NOM,IIO,JJO,KKO
TYPE (MESH_TYPE), POINTER :: M2=>NULL()

OTHER_MESH_LOOP: DO NOM=1,NMESHES
   IF (EVACUATION_ONLY(NOM)) CYCLE OTHER_MESH_LOOP
   M2=>MESHES(NOM)
   IF (XX>=M2%XS .AND. XX<=M2%XF .AND.  YY>=M2%YS .AND. YY<=M2%YF .AND. ZZ>=M2%ZS .AND. ZZ<=M2%ZF) THEN
      XI  = MAX( 1._EB , MIN( REAL(M2%IBAR,EB)+ALMOST_ONE , M2%CELLSI(FLOOR((XX-M2%XS)*M2%RDXINT)) + 1._EB ) )
      YJ  = MAX( 1._EB , MIN( REAL(M2%JBAR,EB)+ALMOST_ONE , M2%CELLSJ(FLOOR((YY-M2%YS)*M2%RDYINT)) + 1._EB ) )
      ZK  = MAX( 1._EB , MIN( REAL(M2%KBAR,EB)+ALMOST_ONE , M2%CELLSK(FLOOR((ZZ-M2%ZS)*M2%RDZINT)) + 1._EB ) )
      IIO = FLOOR(XI)
      JJO = FLOOR(YJ)
      KKO = FLOOR(ZK)
      RETURN
   ENDIF
ENDDO OTHER_MESH_LOOP

NOM = 0

END SUBROUTINE SEARCH_OTHER_MESHES


SUBROUTINE BLOCK_CELL(NM,I1,I2,J1,J2,K1,K2,IVAL,OBST_INDEX)

! Indicate which cells are blocked off

INTEGER :: NM,I1,I2,J1,J2,K1,K2,IVAL,I,J,K,OBST_INDEX,IC
TYPE (MESH_TYPE), POINTER :: M=>NULL()

M => MESHES(NM)
DO K=K1,K2
   DO J=J1,J2
      DO I=I1,I2
         IC = M%CELL_INDEX(I,J,K)
         SELECT CASE(IVAL)
            CASE(0)
               M%SOLID(IC)        = .FALSE.
               M%OBST_INDEX_C(IC) = 0
            CASE(1)
               M%SOLID(IC)        = .TRUE.
               M%OBST_INDEX_C(IC) = OBST_INDEX
         END SELECT
         IF (OBST_INDEX==0) M%EXTERIOR(IC) = .TRUE.
      ENDDO
   ENDDO
ENDDO

END SUBROUTINE BLOCK_CELL


SUBROUTINE ASSIGN_HT3D_WALL_INDICES(NM)

! For each cell with CELL_INDEX=IC in an HT3D solid obstruction, fill in the array MESHES(NM)%WALL_INDEX_HT3D(IC,-3:3) such that:
! WALL_INDEX_HT3D(IC,-3) is the wall index of the bottom of the solid obstruction
! WALL_INDEX_HT3D(IC,-2) is the wall index of the back   of the solid obstruction
! WALL_INDEX_HT3D(IC,-1) is the wall index of the left   of the solid obstruction
! WALL_INDEX_HT3D(IC, 1) is the wall index of the right  of the solid obstruction
! WALL_INDEX_HT3D(IC, 2) is the wall index of the front  of the solid obstruction
! WALL_INDEX_HT3D(IC, 3) is the wall index of the top    of the solid obstruction
! WALL_INDEX_HT3D(IC, 0) is the wall index of the nearest surface of the solid obstruction

INTEGER, INTENT(IN) :: NM
INTEGER :: I,J,K,N,IC,II,JJ,KK,ICN,CELL_COUNT(-3:3)
TYPE(MESH_TYPE), POINTER :: M
TYPE(OBSTRUCTION_TYPE), POINTER :: OB

M => MESHES(NM)

OBST_LOOP: DO N=1,M%N_OBST
   OB => M%OBSTRUCTION(N)
   HT3D_IF: IF (OB%HT3D) THEN
      K_LOOP: DO K=OB%K1+1,OB%K2
         J_LOOP: DO J=OB%J1+1,OB%J2
            I_LOOP: DO I=OB%I1+1,OB%I2
               IC = M%CELL_INDEX(I,J,K)
               IF (.NOT.M%SOLID(IC)) CYCLE I_LOOP
               M%WALL_INDEX_HT3D(IC,:) = 0
               CELL_COUNT = 0
               CELL_COUNT(0) = 1000000

               MARCH_RIGHT: DO II=I+1,M%IBAR
                  ICN = M%CELL_INDEX(II,J,K)
                  CELL_COUNT(1) = CELL_COUNT(1) + 1
                  IF (.NOT.M%SOLID(ICN)) THEN
                     M%WALL_INDEX_HT3D(IC,1) = M%WALL_INDEX(ICN,-1)
                     EXIT MARCH_RIGHT
                  ENDIF
                  IF (II==M%IBAR) CELL_COUNT(1) = 1000000
               ENDDO MARCH_RIGHT

               MARCH_LEFT: DO II=I-1,1,-1
                  ICN = M%CELL_INDEX(II,J,K)
                  CELL_COUNT(-1) = CELL_COUNT(-1) + 1
                  IF (.NOT.M%SOLID(ICN)) THEN
                     M%WALL_INDEX_HT3D(IC,-1) = M%WALL_INDEX(ICN,1)
                     EXIT MARCH_LEFT
                  ENDIF
                  IF (II==1) CELL_COUNT(-1) = 1000000
               ENDDO MARCH_LEFT

               MARCH_FORWARD: DO JJ=J+1,M%JBAR
                  ICN = M%CELL_INDEX(I,JJ,K)
                  CELL_COUNT(2) = CELL_COUNT(2) + 1
                  IF (.NOT.M%SOLID(ICN)) THEN
                     M%WALL_INDEX_HT3D(IC,2) = M%WALL_INDEX(ICN,-2)
                     EXIT MARCH_FORWARD
                  ENDIF
                  IF (JJ==M%JBAR) CELL_COUNT(2) = 1000000
               ENDDO MARCH_FORWARD

               MARCH_BACK: DO JJ=J-1,1,-1
                  ICN = M%CELL_INDEX(I,JJ,K)
                  CELL_COUNT(-2) = CELL_COUNT(-2) + 1
                  IF (.NOT.M%SOLID(ICN)) THEN
                     M%WALL_INDEX_HT3D(IC,-2) = M%WALL_INDEX(ICN,2)
                     EXIT MARCH_BACK
                  ENDIF
                  IF (JJ==1) CELL_COUNT(-2) = 1000000
               ENDDO MARCH_BACK

               MARCH_UP: DO KK=K+1,M%KBAR
                  ICN = M%CELL_INDEX(I,J,KK)
                  CELL_COUNT(3) = CELL_COUNT(3) + 1
                  IF (.NOT.M%SOLID(ICN)) THEN
                     M%WALL_INDEX_HT3D(IC,3) = M%WALL_INDEX(ICN,-3)
                     EXIT MARCH_UP
                  ENDIF
                  IF (KK==M%KBAR) CELL_COUNT(3) = 1000000
               ENDDO MARCH_UP

               MARCH_DOWN: DO KK=K-1,1,-1
                  ICN = M%CELL_INDEX(I,J,KK)
                  CELL_COUNT(-3) = CELL_COUNT(-3) + 1
                  IF (.NOT.M%SOLID(ICN)) THEN
                     M%WALL_INDEX_HT3D(IC,-3) = M%WALL_INDEX(ICN,3)
                     EXIT MARCH_DOWN
                  ENDIF
                  IF (KK==1) CELL_COUNT(-3) = 1000000
               ENDDO MARCH_DOWN

               M%WALL_INDEX_HT3D(IC,0) = M%WALL_INDEX_HT3D(IC,MINLOC(CELL_COUNT,DIM=1)-4)

            ENDDO I_LOOP
         ENDDO J_LOOP
      ENDDO K_LOOP
   ENDIF HT3D_IF
ENDDO OBST_LOOP

END SUBROUTINE ASSIGN_HT3D_WALL_INDICES


LOGICAL FUNCTION INTERIOR(XX,YY,ZZ)

INTEGER NM
REAL(EB), INTENT(IN) :: XX,YY,ZZ

INTERIOR = .FALSE.

DO NM=1,NMESHES
   IF (XX>MESHES(NM)%XS .AND. XX<MESHES(NM)%XF .AND. &
       YY>MESHES(NM)%YS .AND. YY<MESHES(NM)%YF .AND. &
       ZZ>MESHES(NM)%ZS .AND. ZZ<MESHES(NM)%ZF) INTERIOR = .TRUE.
ENDDO

END FUNCTION INTERIOR


SUBROUTINE ASSIGN_PRESSURE_ZONE(NM,XX,YY,ZZ,I_ZONE,I_ZONE_OVERLAP)

! Given the point (XX,YY,ZZ) within Mesh NM, determine all mesh cells within the same pressure zone

REAL(EB), INTENT(IN) :: XX,YY,ZZ
REAL(EB) :: XI,YJ,ZK
INTEGER, INTENT(IN) :: NM,I_ZONE
INTEGER, INTENT(OUT) :: I_ZONE_OVERLAP
INTEGER :: NN,IOR,IC,II,JJ,KK,III,JJJ,KKK,Q_N,IIO,JJO,KKO,NOM
INTEGER, ALLOCATABLE, DIMENSION(:) :: Q_I,Q_J,Q_K
TYPE (MESH_TYPE), POINTER :: M=>NULL()
TYPE (OBSTRUCTION_TYPE), POINTER :: OB=>NULL()

M=>MESHES(NM)
I_ZONE_OVERLAP = 0

IF (XX<M%XS-TWO_EPSILON_EB .OR. XX>M%XF+TWO_EPSILON_EB .OR. &
    YY<M%YS-TWO_EPSILON_EB .OR. YY>M%YF+TWO_EPSILON_EB .OR. &
    ZZ<M%ZS-TWO_EPSILON_EB .OR. ZZ>M%ZF+TWO_EPSILON_EB) RETURN

ALLOCATE(Q_I(M%IBAR*M%JBAR*M%KBAR))
ALLOCATE(Q_J(M%IBAR*M%JBAR*M%KBAR))
ALLOCATE(Q_K(M%IBAR*M%JBAR*M%KBAR))

! Find the cell indices corresponding to the given point

XI  = MAX( 1._EB , MIN( REAL(M%IBAR,EB)+ALMOST_ONE , M%CELLSI(NINT((XX-M%XS)*M%RDXINT)) + 1._EB ) )
YJ  = MAX( 1._EB , MIN( REAL(M%JBAR,EB)+ALMOST_ONE , M%CELLSJ(NINT((YY-M%YS)*M%RDYINT)) + 1._EB ) )
ZK  = MAX( 1._EB , MIN( REAL(M%KBAR,EB)+ALMOST_ONE , M%CELLSK(NINT((ZZ-M%ZS)*M%RDZINT)) + 1._EB ) )
II  = FLOOR(XI)
JJ  = FLOOR(YJ)
KK  = FLOOR(ZK)

! Add the first entry to "queue" of cells that need a pressure zone number

Q_I(1) = II
Q_J(1) = JJ
Q_K(1) = KK
Q_N    = 1
M%PRESSURE_ZONE(II,JJ,KK) = I_ZONE

! Look to all cells adjacent to the starting cell and determine if they are in the ZONE as well.
! Repeat process until all cells are found.

SORT_QUEUE: DO

   IF (Q_N<1) EXIT SORT_QUEUE

   III = Q_I(Q_N)
   JJJ = Q_J(Q_N)
   KKK = Q_K(Q_N)
   IC  = M%CELL_INDEX(III,JJJ,KKK)
   Q_N = Q_N - 1

   SEARCH_LOOP: DO IOR=-3,3

      IF (IOR==0) CYCLE SEARCH_LOOP

      SELECT CASE(IOR)
         CASE(-1)
            IF (III==0) CYCLE SEARCH_LOOP
            II = III-1
            JJ = JJJ
            KK = KKK
         CASE( 1)
            IF (III==M%IBP1) CYCLE SEARCH_LOOP
            II = III+1
            JJ = JJJ
            KK = KKK
         CASE(-2)
            IF (JJJ==0) CYCLE SEARCH_LOOP
            II = III
            JJ = JJJ-1
            KK = KKK
         CASE( 2)
            IF (JJJ==M%JBP1) CYCLE SEARCH_LOOP
            II = III
            JJ = JJJ+1
            KK = KKK
         CASE(-3)
            IF (KKK==0) CYCLE SEARCH_LOOP
            II = III
            JJ = JJJ
            KK = KKK-1
         CASE( 3)
            IF (KKK==M%KBP1) CYCLE SEARCH_LOOP
            II = III
            JJ = JJJ
            KK = KKK+1
      END SELECT

      ! If the cell is outside the computational domain, check if it is in another mesh

      IF (II<1 .OR. II>M%IBAR .OR. JJ<1 .OR. JJ>M%JBAR .OR. KK<1 .OR. KK>M%KBAR) THEN
         CALL SEARCH_OTHER_MESHES(M%XC(II),M%YC(JJ),M%ZC(KK),NOM,IIO,JJO,KKO)
         IF (NOM>0) M%PRESSURE_ZONE(II,JJ,KK) = I_ZONE
         CYCLE SEARCH_LOOP
      ENDIF

      ! Look for thin obstructions bordering the current cell

      DO NN=1,M%N_OBST
         OB=>M%OBSTRUCTION(NN)
         SELECT CASE(IOR)
            CASE(-1)
               IF (II==  OB%I1 .AND. II==  OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
            CASE( 1)
               IF (II-1==OB%I1 .AND. II-1==OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
            CASE(-2)
               IF (JJ==  OB%J1 .AND. JJ==  OB%J2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
            CASE( 2)
               IF (JJ-1==OB%J1 .AND. JJ-1==OB%J2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. KK>OB%K1 .AND. KK<=OB%K2) CYCLE SEARCH_LOOP
            CASE(-3)
               IF (KK==  OB%K1 .AND. KK==  OB%K2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2) CYCLE SEARCH_LOOP
            CASE( 3)
               IF (KK-1==OB%K1 .AND. KK-1==OB%K2 .AND. II>OB%I1 .AND. II<=OB%I2 .AND. JJ>OB%J1 .AND. JJ<=OB%J2) CYCLE SEARCH_LOOP
         END SELECT
      ENDDO

      ! If an obstruction is found, assign its cells the current ZONE, just in case the obstruction is removed

      IC = M%CELL_INDEX(II,JJ,KK)
      IF (M%SOLID(IC) .AND. M%OBST_INDEX_C(IC)>0) THEN
         OB => M%OBSTRUCTION(M%OBST_INDEX_C(IC))
         M%PRESSURE_ZONE(OB%I1+1:OB%I2,OB%J1+1:OB%J2,OB%K1+1:OB%K2) = I_ZONE
         CYCLE SEARCH_LOOP
      ENDIF

      ! If the neighboring cell is not solid, assign the pressure zone

      IF (.NOT.M%SOLID(IC) .AND. M%PRESSURE_ZONE(II,JJ,KK)>0 .AND.  M%PRESSURE_ZONE(II,JJ,KK)/=I_ZONE) THEN
         I_ZONE_OVERLAP = M%PRESSURE_ZONE(II,JJ,KK)
         RETURN
      ENDIF

      IF (.NOT.M%SOLID(IC) .AND. M%PRESSURE_ZONE(II,JJ,KK)<1) THEN
         Q_N      = Q_N+1
         Q_I(Q_N) = II
         Q_J(Q_N) = JJ
         Q_K(Q_N) = KK
         M%PRESSURE_ZONE(II,JJ,KK) = I_ZONE
      ENDIF

   ENDDO SEARCH_LOOP

END DO SORT_QUEUE

DEALLOCATE(Q_I)
DEALLOCATE(Q_J)
DEALLOCATE(Q_K)

END SUBROUTINE ASSIGN_PRESSURE_ZONE


SUBROUTINE GET_N_LAYER_CELLS(DIFFUSIVITY,LAYER_THICKNESS,STRETCH_FACTOR,CELL_SIZE_FACTOR,N_LAYER_CELLS_MAX,N_CELLS,DX_MIN)

! Get number of wall cells in the layer

INTEGER, INTENT(OUT)  :: N_CELLS
INTEGER, INTENT(IN) :: N_LAYER_CELLS_MAX
REAL(EB), INTENT(OUT) :: DX_MIN
REAL(EB), INTENT(IN)  :: DIFFUSIVITY,LAYER_THICKNESS,STRETCH_FACTOR,CELL_SIZE_FACTOR
REAL(EB) :: DDSUM
INTEGER  :: N, I

DX_MIN = 0._EB
IF (LAYER_THICKNESS<=TWO_EPSILON_EB) THEN
   N_CELLS = 0
   DX_MIN = 0._EB
   RETURN
ENDIF

SHRINK_LOOP: DO N=1,N_LAYER_CELLS_MAX-1
   DDSUM = 0._EB
   SUM_LOOP: DO I=1,N
      DDSUM = DDSUM + STRETCH_FACTOR**(MIN(I-1,N-I))
   ENDDO SUM_LOOP
   IF ((LAYER_THICKNESS/DDSUM < CELL_SIZE_FACTOR*SQRT(DIFFUSIVITY)) .OR. (N==N_LAYER_CELLS_MAX-1)) THEN
      N_CELLS = N
      DX_MIN = LAYER_THICKNESS/DDSUM
      EXIT SHRINK_LOOP
   ENDIF
ENDDO SHRINK_LOOP

END SUBROUTINE GET_N_LAYER_CELLS


SUBROUTINE GET_WALL_NODE_COORDINATES(N_CELLS,N_LAYERS,N_LAYER_CELLS,SMALLEST_CELL_SIZE,STRETCH_FACTOR,X_S)

! Get the wall internal coordinates

INTEGER, INTENT(IN)     :: N_CELLS,N_LAYERS, N_LAYER_CELLS(N_LAYERS)
REAL(EB), INTENT(IN)    :: SMALLEST_CELL_SIZE(N_LAYERS),STRETCH_FACTOR(N_LAYERS)
REAL(EB), INTENT(OUT)   :: X_S(0:N_CELLS)

INTEGER I, II, NL
REAL(EB) DX_S

II = 0
X_S(0) = 0._EB
DO NL=1,N_LAYERS
   DO I=1,N_LAYER_CELLS(NL)
      II = II+1
      DX_S = SMALLEST_CELL_SIZE(NL)*STRETCH_FACTOR(NL)**(MIN(I-1,N_LAYER_CELLS(NL)-I))
      X_S(II) = X_S(II-1) + DX_S
   ENDDO
ENDDO

END SUBROUTINE GET_WALL_NODE_COORDINATES



SUBROUTINE GET_WALL_NODE_WEIGHTS(N_CELLS,N_LAYERS,N_LAYER_CELLS, &
         LAYER_THICKNESS,GEOMETRY,X_S,LAYER_DIVIDE,DX,RDX,RDXN,DX_WGT,DXF,DXB,LAYER_INDEX,MF_FRAC,INNER_RADIUS)

! Get the wall internal coordinates

INTEGER, INTENT(IN)     :: N_CELLS, N_LAYERS, N_LAYER_CELLS(N_LAYERS),GEOMETRY
REAL(EB), INTENT(IN)    :: X_S(0:N_CELLS),LAYER_THICKNESS(1:N_LAYERS),LAYER_DIVIDE,INNER_RADIUS
INTEGER, INTENT(OUT)    :: LAYER_INDEX(0:N_CELLS+1)
REAL(EB), INTENT(OUT)   :: DX(1:N_CELLS),RDX(0:N_CELLS+1),RDXN(0:N_CELLS),DX_WGT(0:N_CELLS),DXF,DXB, &
                           MF_FRAC(1:N_CELLS)

INTEGER :: I, II, NL, I_GRAD=0
REAL(EB) :: R, THICKNESS, X_DIVIDE

   THICKNESS = SUM(LAYER_THICKNESS)

   SELECT CASE(GEOMETRY)
   CASE(SURF_CARTESIAN)
      I_GRAD = 0
   CASE(SURF_CYLINDRICAL)
      I_GRAD = 1
   CASE(SURF_SPHERICAL)
      I_GRAD = 2
   END SELECT

   II = 0
   DO NL=1,N_LAYERS
      DO I=1,N_LAYER_CELLS(NL)
         II = II + 1
         LAYER_INDEX(II) = NL
      ENDDO
   ENDDO
   LAYER_INDEX(0) = 1
   LAYER_INDEX(N_CELLS+1) = N_LAYERS
   DXF = X_S(1)       - X_S(0)
   DXB = X_S(N_CELLS) - X_S(N_CELLS-1)

   ! Compute dx_weight for each node (dx_weight is the ratio of cell size to the combined size of the current and next cell)

   DO I=1,N_CELLS-1
      DX_WGT(I) = (X_S(I)-X_S(I-1))/(X_S(I+1)-X_S(I-1))
   ENDDO
   DX_WGT(0)       = 0.5_EB
   DX_WGT(N_CELLS) = 0.5_EB

   ! Compute dx and 1/dx for each node (dx is the distance from cell boundary to cell boundary)

   DO I=1,N_CELLS
      DX(I)  = X_S(I)-X_S(I-1)
      RDX(I) = 1._EB/DX(I)
   ENDDO

   ! Adjust 1/dx_n to 1/rdr for cylindrical case and 1/r^2dr for spaherical

   SELECT CASE(GEOMETRY)
      CASE(SURF_CYLINDRICAL)
         DO I=1,N_CELLS
            R = 0.5_EB*RDX(I)*((INNER_RADIUS+THICKNESS-X_S(I-1))**2-(INNER_RADIUS+THICKNESS-X_S(I))**2)
            RDX(I) = RDX(I)/R**I_GRAD
         ENDDO
      CASE(SURF_SPHERICAL)
         DO I=1,N_CELLS
            R = SQRT(ONTH*RDX(I)*((INNER_RADIUS+THICKNESS-X_S(I-1))**3-(INNER_RADIUS+THICKNESS-X_S(I))**3))
            RDX(I) = RDX(I)/R**I_GRAD
         ENDDO
   END SELECT
   RDX(0)         = RDX(1)
   RDX(N_CELLS+1) = RDX(N_CELLS)

   ! Compute 1/dx_n for each node (dx_n is the distance from cell center to cell center)

   DO I=1,N_CELLS-1
      RDXN(I) = 2._EB/(X_S(I+1)-X_S(I-1))
   ENDDO
   RDXN(0)       = 1._EB/(X_S(1)-X_S(0))
   RDXN(N_CELLS) = 1._EB/(X_S(N_CELLS)-X_S(N_CELLS-1))

   ! Adjust 1/dx_n to r/dr for cylindrical case and r^2/dr for spaherical

   IF (GEOMETRY /= SURF_CARTESIAN) THEN
      DO I=0,N_CELLS
         R = INNER_RADIUS+THICKNESS-X_S(I)
         RDXN(I) = RDXN(I)*R**I_GRAD
      ENDDO
   ENDIF

   ! Compute mass flux fraction array (array numbers indicate the fraction of mass flux that is added to the front

   IF (LAYER_DIVIDE >= REAL(N_LAYERS,EB)) THEN
      MF_FRAC = 1.0_EB
   ELSE
      MF_FRAC = 0._EB

      X_DIVIDE = 0._EB
      DO NL=1,N_LAYERS
         IF (LAYER_DIVIDE>=REAL(NL,EB)) THEN
            X_DIVIDE  = X_DIVIDE + LAYER_THICKNESS(NL)
         ELSE
            X_DIVIDE  = X_DIVIDE + MOD(LAYER_DIVIDE,1.0_EB)*LAYER_THICKNESS(NL)
            EXIT
         ENDIF
      ENDDO

      II = 0
      DIVILOOP: DO NL=1,N_LAYERS
         DO I=1,N_LAYER_CELLS(NL)
            II = II + 1
            IF (X_S(II) < X_DIVIDE) THEN
               MF_FRAC(II) = 1._EB
            ELSEIF (X_S(II-1) < X_DIVIDE) THEN
               MF_FRAC(II) = (X_DIVIDE-X_S(II-1))/DX(II)
               EXIT DIVILOOP
            ENDIF
         ENDDO
      ENDDO DIVILOOP
   ENDIF

END SUBROUTINE GET_WALL_NODE_WEIGHTS



SUBROUTINE GET_INTERPOLATION_WEIGHTS(N_LAYERS,NWP,NWP_NEW,N_LAYER_CELLS,N_LAYER_CELLS_NEW,X_S,X_S_NEW,INT_WGT)

INTEGER, INTENT(IN)  :: N_LAYERS,NWP,NWP_NEW,N_LAYER_CELLS(N_LAYERS),N_LAYER_CELLS_NEW(N_LAYERS)
REAL(EB), INTENT(IN) :: X_S(0:NWP), X_S_NEW(0:NWP_NEW)
REAL(EB), INTENT(OUT) :: INT_WGT(:,:)

REAL(EB) XUP,XLOW,XUP_NEW,XLOW_NEW,DX_NEW
INTEGER I, J, II, JJ, I_BASE, J_BASE, J_OLD,N
II = 0
JJ = 0
I_BASE = 0
J_BASE = 0

INT_WGT = 0._EB
DO N = 1,N_LAYERS
   J_OLD = 1
   DO I = 1,N_LAYER_CELLS_NEW(N)
      II       = I_BASE + I
      XUP_NEW  = X_S_NEW(II)
      XLOW_NEW = X_S_NEW(II-1)
      DX_NEW   = XUP_NEW - XLOW_NEW
   !  XUP_NEW  = X_S_NEW(NWP_NEW)-X_S_NEW(II-1)
   !  XLOW_NEW = X_S_NEW(NWP_NEW)-X_S_NEW(II)
   !  DX_NEW   = XUP_NEW**3 - XLOW_NEW**3
      DO J = J_OLD,N_LAYER_CELLS(N)
         JJ = J_BASE + J
         XUP =  X_S(JJ)
         XLOW = X_S(JJ-1)
   !     XUP =  X_S(NWP)-X_S(JJ-1)
   !     XLOW = X_S(NWP)-X_S(JJ)
         INT_WGT(II,JJ) = (MIN(XUP,XUP_NEW)-MAX(XLOW,XLOW_NEW))/DX_NEW
   !     INT_WGT(II,JJ) = (MIN(XUP,XUP_NEW)**3-MAX(XLOW,XLOW_NEW)**3)/DX_NEW
         IF (XUP >= XUP_NEW) EXIT
      ENDDO
      J_OLD = J
   ENDDO
   I_BASE = I_BASE + N_LAYER_CELLS_NEW(N)
   J_BASE = J_BASE + N_LAYER_CELLS(N)
ENDDO

END SUBROUTINE GET_INTERPOLATION_WEIGHTS



SUBROUTINE INTERPOLATE_WALL_ARRAY(N_CELLS,NWP,NWP_NEW,INT_WGT,ARR)

INTEGER, INTENT(IN)  :: N_CELLS,NWP,NWP_NEW
REAL(EB), INTENT(IN) :: INT_WGT(:,:)
REAL(EB) ARR(N_CELLS),TMP(N_CELLS)

INTEGER I,J

TMP = ARR
ARR = 0._EB
DO I = 1,NWP_NEW
DO J = 1,NWP
   ARR(I) = ARR(I) + INT_WGT(I,J)*TMP(J)
ENDDO
ENDDO

END SUBROUTINE INTERPOLATE_WALL_ARRAY


SUBROUTINE RANDOM_RECTANGLE(XX,YY,ZZ,X1,X2,Y1,Y2,Z1,Z2)

! Choose a random point (XX,YY,ZZ) from within a rectangular volume bounded by X1, X2, ...

REAL(EB), INTENT(IN) :: X1,X2,Y1,Y2,Z1,Z2
REAL(EB), INTENT(OUT) :: XX,YY,ZZ
REAL(EB) :: RN

CALL RANDOM_NUMBER(RN)
XX = X1 + RN*(X2-X1)
CALL RANDOM_NUMBER(RN)
YY = Y1 + RN*(Y2-Y1)
CALL RANDOM_NUMBER(RN)
ZZ = Z1 + RN*(Z2-Z1)

END SUBROUTINE RANDOM_RECTANGLE


SUBROUTINE RANDOM_CONE(NM,XX,YY,ZZ,X0,Y0,Z0,RR0,HH0)

! Choose a random point (XX,YY,ZZ) from within a vertically oriented cone

INTEGER, INTENT(IN) :: NM
REAL(EB), INTENT(IN) :: X0,Y0,Z0,RR0,HH0
REAL(EB), INTENT(OUT) :: XX,YY,ZZ
REAL(EB) :: THETA,RR,RN
TYPE (MESH_TYPE), POINTER :: M

M => MESHES(NM)
SEARCH_LOOP: DO
   CALL RANDOM_NUMBER(RN)
   ZZ = Z0 + HH0*RN
   CALL RANDOM_NUMBER(RN)
   THETA = TWOPI*RN
   CALL RANDOM_NUMBER(RN)
   RR = SQRT(RN)*RR0
   XX = X0 + RR*COS(THETA)
   YY = Y0 + RR*SIN(THETA)
   IF (RR**2>(RR0*(1._EB-(ZZ-Z0)/HH0))**2) CYCLE SEARCH_LOOP
   IF (XX>=M%XS .AND. XX<=M%XF .AND. YY>=M%YS .AND. YY<=M%YF .AND. ZZ>=M%ZS .AND. ZZ<=M%ZF) EXIT SEARCH_LOOP
ENDDO SEARCH_LOOP

END SUBROUTINE RANDOM_CONE


SUBROUTINE RANDOM_RING(NM,XX,YY,X0,Y0,RR0)

! Choose a random point (XX,YY) on a horizontally-oriented ring

INTEGER, INTENT(IN) :: NM
REAL(EB), INTENT(IN) :: X0,Y0,RR0
REAL(EB), INTENT(OUT) :: XX,YY
REAL(EB) :: THETA,RN
TYPE (MESH_TYPE), POINTER :: M

M => MESHES(NM)
SEARCH_LOOP: DO
   CALL RANDOM_NUMBER(RN)
   THETA = TWOPI*RN
   XX = X0 + RR0*COS(THETA)
   YY = Y0 + RR0*SIN(THETA)
   IF (XX>=M%XS .AND. XX<=M%XF .AND. YY>=M%YS .AND. YY<=M%YF) EXIT SEARCH_LOOP
ENDDO SEARCH_LOOP

END SUBROUTINE RANDOM_RING


SUBROUTINE UNIFORM_RING(NM,XX,YY,X0,Y0,RR0,NP,NP_TOTAL)

! Choose a point (XX,YY) uniformly on a horizontally-oriented ring

INTEGER, INTENT(IN) :: NM,NP,NP_TOTAL
REAL(EB), INTENT(IN) :: X0,Y0,RR0
REAL(EB), INTENT(OUT) :: XX,YY
REAL(EB) :: THETA
TYPE (MESH_TYPE), POINTER :: M

M => MESHES(NM)
SEARCH_LOOP: DO
   THETA = TWOPI*REAL(NP,EB)/REAL(NP_TOTAL,EB)
   XX = X0 + RR0*COS(THETA)
   YY = Y0 + RR0*SIN(THETA)
   IF (XX>=M%XS .AND. XX<=M%XF .AND. YY>=M%YS .AND. YY<=M%YF) EXIT SEARCH_LOOP
ENDDO SEARCH_LOOP

END SUBROUTINE UNIFORM_RING


REAL(EB) FUNCTION CIRCLE_CELL_INTERSECTION_AREA(X0,Y0,RAD,X1,X2,Y1,Y2)

! Estimate area of intersection of circle with origin (X0,Y0), radius, RAD, and rectangle (X1,Y1,X2,Y2)

REAL(EB), INTENT(IN) :: X0,Y0,RAD,X1,X2,Y1,Y2
INTEGER :: NN,II,JJ
REAL(EB) :: DELTA_AREA,XX,YY

CIRCLE_CELL_INTERSECTION_AREA = 0._EB
NN = 50
DELTA_AREA = (X2-X1)*(Y2-Y1)/REAL(NN,EB)**2

DO JJ=1,NN
   DO II=1,NN
      XX = X1 + (X2-X1)*(II-0.5_EB)/REAL(NN,EB)
      YY = Y1 + (Y2-Y1)*(JJ-0.5_EB)/REAL(NN,EB)
      IF ( ((XX-X0)**2+(YY-Y0)**2)<RAD**2 ) CIRCLE_CELL_INTERSECTION_AREA = CIRCLE_CELL_INTERSECTION_AREA + DELTA_AREA
   ENDDO
ENDDO

END FUNCTION CIRCLE_CELL_INTERSECTION_AREA


REAL(EB) FUNCTION CONE_MESH_INTERSECTION_VOLUME(NM,X0,Y0,Z0,RR0,HH0)

! Calculate volume of the interaction of MESH NM with a cone

INTEGER, INTENT(IN) :: NM
REAL(EB), INTENT(IN) :: X0,Y0,Z0,RR0,HH0
INTEGER :: I,J,K,NX,NY,NZ
REAL(EB) :: XX,YY,ZZ,R2,X_MIN,X_MAX,Y_MIN,Y_MAX,Z_MIN,Z_MAX,DX,DY,DZ
TYPE (MESH_TYPE), POINTER :: M

CONE_MESH_INTERSECTION_VOLUME = 0._EB
M => MESHES(NM)

X_MIN = MAX(M%XS,X0-RR0)
X_MAX = MIN(M%XF,X0+RR0)
Y_MIN = MAX(M%YS,Y0-RR0)
Y_MAX = MIN(M%YF,Y0+RR0)
Z_MIN = MAX(M%ZS,Z0)
Z_MAX = MIN(M%ZF,Z0+HH0)
IF (X_MAX<=X_MIN .OR. Y_MAX<=Y_MIN .OR. Z_MAX<=Z_MIN) RETURN
NX = CEILING(10*(X_MAX-X_MIN))
NY = CEILING(10*(Y_MAX-Y_MIN))
NZ = CEILING(10*(Z_MAX-Z_MIN))
DX = (X_MAX-X_MIN)/REAL(NX,EB)
DY = (Y_MAX-Y_MIN)/REAL(NY,EB)
DZ = (Z_MAX-Z_MIN)/REAL(NZ,EB)

DO K=1,NZ
   ZZ = Z_MIN + (K-0.5_EB)*DZ
   IF (ZZ<Z0 .OR. ZZ>Z0+HH0) CYCLE
   DO J=1,NY
      YY = Y_MIN + (J-0.5_EB)*DY
      DO I=1,NX
         XX = X_MIN + (I-0.5_EB)*DX
         R2 = (XX-X0)**2+(YY-Y0)**2
         IF (R2<(RR0*(1._EB-(ZZ-Z0)/HH0))**2) &
            CONE_MESH_INTERSECTION_VOLUME = CONE_MESH_INTERSECTION_VOLUME + DX*DY*DZ
      ENDDO
   ENDDO
ENDDO

END FUNCTION CONE_MESH_INTERSECTION_VOLUME


SUBROUTINE TRANSFORM_COORDINATES(X,Y,Z,MOVE_INDEX)

! Transform (X,Y,Z) using the parameters of MOVE(MOVE_INDEX)

INTEGER, INTENT(IN) :: MOVE_INDEX
REAL(EB), INTENT(INOUT) :: X,Y,Z
REAL(EB) :: M(3,3),UP(3,1),S(3,3),UUT(3,3),IDENTITY(3,3),X_VECTOR(3,1),X_VECTOR_0(3,1)
TYPE(MOVEMENT_TYPE), POINTER :: MV

MV => MOVEMENT(MOVE_INDEX)

X_VECTOR = RESHAPE( (/X,Y,Z/),(/3,1/) )
X_VECTOR_0 = RESHAPE( (/MV%X0,MV%Y0,MV%Z0/),(/3,1/) )
UP = RESHAPE(MV%AXIS,(/3,1/))
S =  RESHAPE( (/  0.0_EB, -UP(3,1),  UP(2,1),&
                 UP(3,1),   0.0_EB, -UP(1,1),&
                -UP(2,1),  UP(1,1),  0.0_EB  /),(/3,3/))
UUT = MATMUL(UP,TRANSPOSE(UP))
IDENTITY = RESHAPE ((/ 1.0_EB,0.0_EB,0.0_EB,&
                       0.0_EB,1.0_EB,0.0_EB,&
                       0.0_EB,0.0_EB,1.0_EB /),(/3,3/))
M = UUT + COS(MV%ROTATION_ANGLE*DEG2RAD)*(IDENTITY - UUT) + SIN(MV%ROTATION_ANGLE*DEG2RAD)*S

X_VECTOR = X_VECTOR_0 + MATMUL(M,X_VECTOR-X_VECTOR_0)
X = X_VECTOR(1,1) + MV%DX
Y = X_VECTOR(2,1) + MV%DY
Z = X_VECTOR(3,1) + MV%DZ

END SUBROUTINE TRANSFORM_COORDINATES

END MODULE GEOMETRY_FUNCTIONS


MODULE MATH_FUNCTIONS

USE PRECISION_PARAMETERS
IMPLICIT NONE

CONTAINS

SUBROUTINE UPDATE_HISTOGRAM(NBINS,LIMITS,COUNTS,VAL,WEIGHT)
INTEGER,INTENT(IN)::NBINS
REAL(EB), INTENT(IN)::LIMITS(2),VAL,WEIGHT
REAL(EB), INTENT(INOUT) :: COUNTS(NBINS)
INTEGER::IND=0
IND=MIN(NBINS,MAX(CEILING((VAL-LIMITS(1))/(LIMITS(2)-LIMITS(1))*NBINS),1))
COUNTS(IND)=COUNTS(IND)+WEIGHT
END SUBROUTINE UPDATE_HISTOGRAM


REAL(EB) FUNCTION AFILL(A111,A211,A121,A221,A112,A212,A122,A222,P,R,S)

! Linear interpolation function that returns the weighted average of values at 8 corners of a parallelpiped. A111, etc, are the
! eight values, and P, R, S are fractions of the distance from the origin in the x, y, and z directions.

REAL(EB) :: A111,A211,A121,A221,A112,A212,A122,A222,P,R,S,PP,RR,SS

PP = 1._EB-P
RR = 1._EB-R
SS = 1._EB-S
AFILL = ((PP*A111+P*A211)*RR + (PP*A121+P*A221)*R)*SS + ((PP*A112+P*A212)*RR + (PP*A122+P*A222)*R)*S

END FUNCTION AFILL


REAL(EB) FUNCTION AFILL2(A,I,J,K,P,R,S)

! Linear interpolation function. Same as AFILL, only it reads in entire array.

REAL(EB), INTENT(IN), DIMENSION(0:,0:,0:) :: A
INTEGER, INTENT(IN) :: I,J,K
REAL(EB) A111,A211,A121,A221,A112,A212,A122,A222,P,R,S,PP,RR,SS

A111 = A(I,J,K)
A211 = A(I+1,J,K)
A121 = A(I,J+1,K)
A221 = A(I+1,J+1,K)
A112 = A(I,J,K+1)
A212 = A(I+1,J,K+1)
A122 = A(I,J+1,K+1)
A222 = A(I+1,J+1,K+1)
PP = 1._EB-P
RR = 1._EB-R
SS = 1._EB-S
AFILL2 = ((PP*A111+P*A211)*RR+(PP*A121+P*A221)*R)*SS+ ((PP*A112+P*A212)*RR+(PP*A122+P*A222)*R)*S

END FUNCTION AFILL2


REAL(EB) FUNCTION POLYVAL(N,TEMP,COEF)
! Calculate the value of polynomial function.
INTEGER N,I
REAL(EB) TEMP, COEF(N), VAL
VAL = 0._EB
DO I=1,N
   VAL  = VAL  + COEF(I)*TEMP**(I-1)
ENDDO
POLYVAL = VAL
END FUNCTION POLYVAL


SUBROUTINE GET_RAMP_INDEX(ID,TYPE,RAMP_INDEX,DUPLICATE_RAMP)

USE GLOBAL_CONSTANTS, ONLY: N_RAMP,RAMP_ID,RAMP_TYPE
USE MEMORY_FUNCTIONS, ONLY: REALLOCATE_CHARACTER_ARRAY
CHARACTER(*), INTENT(IN) :: ID,TYPE
INTEGER, INTENT(OUT) :: RAMP_INDEX
INTEGER :: NR
LOGICAL, INTENT(IN), OPTIONAL :: DUPLICATE_RAMP

IF (ID=='null') THEN
   RAMP_INDEX = 0
   RETURN
ENDIF

IF (.NOT.PRESENT(DUPLICATE_RAMP)) THEN
   SEARCH: DO NR=1,N_RAMP
      IF (ID==RAMP_ID(NR)) THEN
         RAMP_INDEX = NR
         RETURN
      ENDIF
   ENDDO SEARCH
ENDIF

IF (N_RAMP>=SIZE(RAMP_ID)) THEN
   RAMP_ID   => REALLOCATE_CHARACTER_ARRAY(RAMP_ID,  LABEL_LENGTH,1,SIZE(RAMP_ID)+100)
   RAMP_TYPE => REALLOCATE_CHARACTER_ARRAY(RAMP_TYPE,LABEL_LENGTH,1,SIZE(RAMP_ID)+100)
ENDIF

N_RAMP                = N_RAMP + 1
RAMP_INDEX            = N_RAMP
RAMP_ID(RAMP_INDEX)   = ID
RAMP_TYPE(RAMP_INDEX) = TYPE

END SUBROUTINE GET_RAMP_INDEX


SUBROUTINE GET_TABLE_INDEX(ID,TYPE,TABLE_INDEX)

USE GLOBAL_CONSTANTS, ONLY: N_TABLE,TABLE_ID,TABLE_TYPE
CHARACTER(*), INTENT(IN) :: ID
INTEGER, INTENT(IN) :: TYPE
INTEGER, INTENT(OUT) :: TABLE_INDEX
INTEGER :: NT

IF (ID=='null') THEN
   TABLE_INDEX = 0
   RETURN
ENDIF

SEARCH: DO NT=1,N_TABLE
   IF (ID==TABLE_ID(NT)) THEN
      TABLE_INDEX = NT
      RETURN
   ENDIF
ENDDO SEARCH

N_TABLE                = N_TABLE + 1
TABLE_INDEX            = N_TABLE
TABLE_ID(TABLE_INDEX)   = ID
TABLE_TYPE(TABLE_INDEX) = TYPE

END SUBROUTINE GET_TABLE_INDEX


REAL(EB) FUNCTION EVALUATE_RAMP(RAMP_INPUT,TAU,RAMP_INDEX)

! General time ramp up

USE TYPES, ONLY: RAMPS
USE DEVICE_VARIABLES, ONLY: DEVICE
USE CONTROL_VARIABLES, ONLY: CONTROL
REAL(EB), INTENT(IN) :: RAMP_INPUT,TAU
REAL(EB):: RAMP_POSITION
INTEGER, INTENT(IN)  :: RAMP_INDEX

SELECT CASE(RAMP_INDEX)
   CASE(-2)
      EVALUATE_RAMP = MAX( TANH(RAMP_INPUT/TAU), 0._EB )
   CASE(-1)
      EVALUATE_RAMP = MIN( (RAMP_INPUT/TAU)**2 , 1.0_EB )
   CASE( 0)
      EVALUATE_RAMP = 1._EB
   CASE(1:)
      IF (RAMPS(RAMP_INDEX)%DEVC_INDEX > 0) THEN
         RAMP_POSITION = &
            MAX(0._EB,MIN(RAMPS(RAMP_INDEX)%SPAN,DEVICE(RAMPS(RAMP_INDEX)%DEVC_INDEX)%SMOOTHED_VALUE - RAMPS(RAMP_INDEX)%T_MIN))
      ELSEIF (RAMPS(RAMP_INDEX)%CTRL_INDEX > 0) THEN
         RAMP_POSITION = &
            MAX(0._EB,MIN(RAMPS(RAMP_INDEX)%SPAN,CONTROL(RAMPS(RAMP_INDEX)%CTRL_INDEX)%INSTANT_VALUE - RAMPS(RAMP_INDEX)%T_MIN))
      ELSE
         RAMP_POSITION = &
            MAX(0._EB,MIN(RAMPS(RAMP_INDEX)%SPAN,RAMP_INPUT - RAMPS(RAMP_INDEX)%T_MIN))
      ENDIF
      EVALUATE_RAMP = RAMPS(RAMP_INDEX)%INTERPOLATED_DATA(NINT(RAMP_POSITION*RAMPS(RAMP_INDEX)%RDT))
END SELECT

END FUNCTION EVALUATE_RAMP


REAL(EB) FUNCTION IERFC(Y)

! Inverse of ERFC(Y)

REAL(EB), INTENT(IN) :: Y
REAL(EB) :: QA,QB,QC,QD,Q0,Q1,Q2,Q3,Q4,PA,PB,P0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10, &
            P11,P12,P13,P14,P15,P16,P17,P18,P19,P20,P21,P22,X,S,T,U,W,Z
PARAMETER ( &
QA = 9.16461398268964-01_EB, &
QB = 2.31729200323405-01_EB, &
QC = 4.88826640273108-01_EB, &
QD = 1.24610454613712-01_EB, &
Q0 = 4.99999303439796-01_EB, &
Q1 = 1.16065025341614-01_EB, &
Q2 = 1.50689047360223-01_EB, &
Q3 = 2.69999308670029-01_EB, &
Q4 = -7.28846765585675-02_EB)
PARAMETER ( &
PA = 3.97886080735226000+00_EB, &
PB = 1.20782237635245222-01_EB, &
P0 = 2.44044510593190935-01_EB, &
P1 = 4.34397492331430115-01_EB, &
P2 = 6.86265948274097816-01_EB, &
P3 = 9.56464974744799006-01_EB, &
P4 = 1.16374581931560831+00_EB, &
P5 = 1.21448730779995237+00_EB, &
P6 = 1.05375024970847138+00_EB, &
P7 = 7.13657635868730364-01_EB, &
P8 = 3.16847638520135944-01_EB, &
P9 = 1.47297938331485121-02_EB, &
P10 = -1.05872177941595488-01_EB, &
P11 = -7.43424357241784861-02_EB)
PARAMETER ( &
P12 = 2.20995927012179067-03_EB, &
P13 = 3.46494207789099922-02_EB, &
P14 = 1.42961988697898018-02_EB, &
P15 = -1.18598117047771104-02_EB, &
P16 = -1.12749169332504870-02_EB, &
P17 = 3.39721910367775861-03_EB, &
P18 = 6.85649426074558612-03_EB, &
P19 = -7.71708358954120939-04_EB, &
P20 = -3.51287146129100025-03_EB, &
P21 = 1.05739299623423047-04_EB, &
P22 = 1.12648096188977922-03_EB)

Z = Y
IF (Y  > 1._EB) Z = 2._EB - Y
W = QA - LOG(Z)
U = SQRT(W)
S = (QC + LOG(U)) / W
T = 1._EB / (U + QB)

X = U * (1._EB - S * (0.5_EB + S * QD)) - ((((Q4 * T + Q3) * T + Q2) * T + Q1) * T + Q0) * T
T = PA / (PA + X)
U = T - 0.5_EB

S = (((((((((P22 * U + P21) * U + P20) * U + P19) * U + P18) * U + P17) * U + P16) * U + P15) * U + P14) * U + P13) * U + P12

S = ((((((((((((S * U + P11) * U + P10) * U +  P9) * U + P8) * U + P7) * U + P6) * U + P5) * U + P4) * U + P3) &
    * U + P2) * U + P1) * U + P0) * T - Z * EXP(X * X - PB)

X = X + S * (1._EB + X * S)

IF (Y > 1._EB) X = -X

IERFC = X

END FUNCTION IERFC


SUBROUTINE GAUSSJ(A,N,NP,B,M,MP,IERROR)

! Solve a linear system of equations with Gauss-Jordon elimination
! Source: Press et al. "Numerical Recipes"

INTEGER :: M,MP,N,NP,I,ICOL=0,IROW=0,J,K,L,LL,INDXC(NP),INDXR(NP),IPIV(NP)
REAL(EB) :: A(NP,NP),B(NP,MP),BIG,DUM,PIVINV
INTEGER, INTENT(OUT) :: IERROR

IERROR = 0
IPIV(1:N) = 0

DO I=1,N
   BIG = 0._EB
   DO J=1,N
      IF (IPIV(J)/=1) THEN
         DO K=1,N
            IF (IPIV(K)==0) THEN
               IF (ABS(A(J,K))>=BIG) THEN
                  BIG = ABS(A(J,K))
                  IROW = J
                  ICOL = K
               ENDIF
            ELSE IF (IPIV(K)>1) THEN
               IERROR = 103   ! Singular matrix in gaussj
               RETURN
            ENDIF
         ENDDO
      ENDIF
   ENDDO
   IPIV(ICOL) = IPIV(ICOL) + 1
   IF (IROW/=ICOL) THEN
      DO L=1,N
         DUM = A(IROW,L)
         A(IROW,L) = A(ICOL,L)
         A(ICOL,L) = DUM
      ENDDO
      DO L=1,M
         DUM = B(IROW,L)
         B(IROW,L) = B(ICOL,L)
         B(ICOL,L) = DUM
      ENDDO
   ENDIF
   INDXR(I) = IROW
   INDXC(I) = ICOL
   IF (ABS(A(ICOL,ICOL))<=TWO_EPSILON_EB) THEN
      IERROR = 103  ! Singular matrix in gaussj
      RETURN
      ENDIF
   PIVINV = 1._EB/A(ICOL,ICOL)
   A(ICOL,ICOL) = 1._EB
   A(ICOL,1:N) = A(ICOL,1:N) * PIVINV
   B(ICOL,1:M) = B(ICOL,1:M) * PIVINV
   DO LL=1,N
      IF (LL/=ICOL) THEN
         DUM = A(LL,ICOL)
         A(LL,ICOL) = 0._EB
         A(LL,1:N) = A(LL,1:N) - A(ICOL,1:N)*DUM
         B(LL,1:M) = B(LL,1:M) - B(ICOL,1:M)*DUM
      ENDIF
   ENDDO
ENDDO
DO L=N,1,-1
   IF (INDXR(L)/=INDXC(L)) THEN
      DO K=1,N
         DUM = A(K,INDXR(L))
         A(K,INDXR(L)) = A(K,INDXC(L))
         A(K,INDXC(L)) = DUM
      ENDDO
   ENDIF
ENDDO

END SUBROUTINE GAUSSJ


SUBROUTINE INTERPOLATE1D(X,Y,XI,ANS)

REAL(EB), INTENT(IN), DIMENSION(:) :: X, Y
REAL(EB), INTENT(IN) :: XI
REAL(EB), INTENT(OUT) :: ANS
INTEGER I, UX,LX

UX = UBOUND(X,1)
LX = LBOUND(X,1)

IF (XI <= X(LX)) THEN
   ANS = Y(LX)
ELSEIF (XI >= X(UX)) THEN
   ANS = Y(UX)
ELSE
   L1: DO I=LX,UX-1
      IF (ABS(XI -X(I)) <= SPACING(X(I))) THEN
         ANS = Y(I)
         EXIT L1
      ELSEIF (X(I+1)>XI) THEN
         ANS = Y(I)+(XI-X(I))/(X(I+1)-X(I)) * (Y(I+1)-Y(I))
         EXIT L1
      ENDIF
   ENDDO L1
ENDIF

END SUBROUTINE INTERPOLATE1D


SUBROUTINE INTERPOLATE1D_UNIFORM(LOWER,X,XI,ANS)

INTEGER, INTENT(IN) :: LOWER
REAL(EB), INTENT(IN), DIMENSION(LOWER:) :: X
REAL(EB), INTENT(IN) :: XI
REAL(EB), INTENT(OUT) :: ANS
INTEGER I, UX,LX
REAL(EB) :: FRAC

UX = UBOUND(X,1)
LX = LBOUND(X,1)

IF (XI <= LX) THEN
   ANS = X(LX)
ELSEIF (XI >= UX) THEN
   ANS = X(UX)
ELSE
   I = INT(XI)
   FRAC = XI - REAL(I,EB)
   ANS = X(I) + FRAC*(X(I+1)-X(I))
ENDIF

END SUBROUTINE INTERPOLATE1D_UNIFORM


SUBROUTINE INTERPOLATE2D(TABLE_INDEX,XI,YI,ANS)
USE TYPES, ONLY: TABLES,TABLES_TYPE
INTEGER, INTENT(IN) :: TABLE_INDEX
REAL(EB), INTENT(IN) :: XI,YI
REAL(EB), INTENT(OUT) :: ANS
REAL(EB) :: XL,XU,FRAC
INTEGER :: I,J
TYPE (TABLES_TYPE), POINTER :: TA

TA => TABLES(TABLE_INDEX)

! Do 1D for X edges of table
IF (XI <= TA%LX) THEN
   CALL INTERPOLATE1D(TA%Y,TA%Z(1,:),YI,ANS)
   RETURN
ELSEIF (XI >= TA%UX) THEN
   CALL INTERPOLATE1D(TA%Y,TA%Z(TA%NUMBER_ROWS,:),YI,ANS)
   RETURN
ENDIF

! Do 1D for Y edges of table
IF (YI <= TA%LY) THEN
   CALL INTERPOLATE1D(TA%X,TA%Z(:,1),XI,ANS)
   RETURN
ELSEIF(YI >= TA%UY) THEN
   CALL INTERPOLATE1D(TA%X,TA%Z(:,TA%NUMBER_COLUMNS),XI,ANS)
   RETURN
ENDIF

! Search for column, do 1D if exact
DO J = 1, TA%NUMBER_COLUMNS
   IF (ABS(YI-TA%Y(J)) < SPACING(TA%Y(J))) THEN
      CALL INTERPOLATE1D(TA%X,TA%Z(:,J),XI,ANS)
      RETURN
   ENDIF
   IF (TA%Y(J) > YI) EXIT
ENDDO

! Search for row, do 1D if exact
DO I = 1, TA%NUMBER_ROWS
   IF(ABS(XI-TA%X(I)) < SPACING(TA%X(I))) THEN
      CALL INTERPOLATE1D(TA%Y,TA%Z(I,:),YI,ANS)
      RETURN
   ENDIF
   IF (TA%X(I) > XI) EXIT
ENDDO

FRAC = (XI-TA%X(I-1))/(TA%X(I)-TA%X(I-1))
XL = FRAC*(TA%Z(I,J-1)-TA%Z(I-1,J-1))+TA%Z(I-1,J-1)
XU = FRAC*(TA%Z(I,J  )-TA%Z(I-1,J  ))+TA%Z(I-1,J)
ANS = (YI-TA%Y(J-1))/(TA%Y(J)-TA%Y(J-1))*(XU-XL)+XL

END SUBROUTINE INTERPOLATE2D


REAL(EB) FUNCTION NORMAL(MEAN,SIGMA)

! Randomly choose a point from a normal distribution

REAL(EB) :: MEAN,SIGMA,TMP,FAC,GSAVE,RSQ,R1,R2
REAL     :: RN
INTEGER :: FLAG
SAVE FLAG,GSAVE
DATA FLAG /0/
IF (FLAG==0) THEN
   RSQ=2.0_EB
   DO WHILE(RSQ>=1.0_EB.OR.RSQ==0.0_EB)
      CALL RANDOM_NUMBER(RN)
      R1=2.0_EB*REAL(RN,EB)-1.0_EB
      CALL RANDOM_NUMBER(RN)
      R2=2.0_EB*REAL(RN,EB)-1.0_EB
      RSQ=R1*R1+R2*R2
   ENDDO
   FAC=SQRT(-2.0_EB*LOG(RSQ)/RSQ)
   GSAVE=R1*FAC
   TMP=R2*FAC
   FLAG=1
ELSE
   TMP=GSAVE
   FLAG=0
ENDIF
NORMAL=TMP*SIGMA+MEAN
RETURN

END FUNCTION NORMAL


SUBROUTINE CROSS_PRODUCT(C,A,B)

! C = A x B

REAL(EB), INTENT(IN) :: A(3),B(3)
REAL(EB), INTENT(OUT) :: C(3)

C(1) = A(2)*B(3)-A(3)*B(2)
C(2) = A(3)*B(1)-A(1)*B(3)
C(3) = A(1)*B(2)-A(2)*B(1)

END SUBROUTINE CROSS_PRODUCT


SUBROUTINE RANDOM_CHOICE(CDF,VAR,NPTS,CHOICE)

! Randomly choose a value from the distribution with given CDF

INTEGER,  INTENT(IN)  :: NPTS
REAL(EB), INTENT(IN)  :: CDF(0:NPTS),VAR(0:NPTS)
REAL(EB), INTENT(OUT) :: CHOICE
INTEGER  :: IT
REAL(EB) :: CFRAC,A,B
REAL(EB) :: RN
REAL     :: RN2

CALL RANDOM_NUMBER(RN2)
RN = REAL(RN2,EB)
A = MINVAL(CDF)
B = MAXVAL(CDF)
RN = A + (B-A)*RN

CDF_LOOP: DO IT=1,NPTS
   IF (CDF(IT) > RN) THEN
      CFRAC  = (RN-CDF(IT-1))/(CDF(IT)-CDF(IT-1))
      CHOICE = VAR(IT-1) + (VAR(IT)-VAR(IT-1))*CFRAC
      EXIT CDF_LOOP
   ENDIF
ENDDO CDF_LOOP

END SUBROUTINE RANDOM_CHOICE


REAL(EB) FUNCTION MINMOD2(X,Y)
REAL(EB), INTENT(IN) :: X,Y
MINMOD2 = 0.5_EB*(SIGN(1._EB,X)+SIGN(1._EB,Y))*MIN(ABS(X),ABS(Y))
END FUNCTION MINMOD2


REAL(EB) FUNCTION MINMOD4(W,X,Y,Z)
REAL(EB), INTENT(IN) :: W,X,Y,Z
MINMOD4 = 0.125_EB*(SIGN(1._EB,W)+SIGN(1._EB,X))* &
          ABS( (SIGN(1._EB,W)+SIGN(1._EB,Y))*(SIGN(1._EB,W)+SIGN(1._EB,Z)) )* &
          MIN(ABS(W),ABS(X),ABS(Y),ABS(Z))
END FUNCTION MINMOD4


SUBROUTINE BOX_MULLER(Z0,Z1)

! Generate pairs of normally distributed pseudo-random numbers with zero mean and unit variance based on the
! Box-Muller transformation.

REAL(EB), INTENT(OUT) :: Z0,Z1
REAL(EB) :: U1,U2,A

CALL RANDOM_NUMBER(U1)
CALL RANDOM_NUMBER(U2)
A = SQRT(-2._EB*LOG(U1))
Z0 = A*COS(TWOPI*U2)
Z1 = A*SIN(TWOPI*U2)

END SUBROUTINE BOX_MULLER


END MODULE MATH_FUNCTIONS


MODULE PHYSICAL_FUNCTIONS

! Functions for physical quantities

USE PRECISION_PARAMETERS
USE GLOBAL_CONSTANTS
USE MESH_VARIABLES
IMPLICIT NONE

CONTAINS


LOGICAL FUNCTION IS_REALIZABLE(ZZ_IN)

REAL(EB), INTENT(IN) :: ZZ_IN(1:N_TRACKED_SPECIES)
REAL(EB), PARAMETER:: ZERO_MINUS=-EPSILON(0._FB),ONE_PLUS=1._EB+EPSILON(1._FB)

IF (ANY(ZZ_IN<ZERO_MINUS) .OR. SUM(ZZ_IN)>ONE_PLUS) THEN
   IS_REALIZABLE=.FALSE.
ELSE
   IS_REALIZABLE=.TRUE.
ENDIF

END FUNCTION IS_REALIZABLE


SUBROUTINE GET_REALIZABLE_MF(ZZ_GET)

REAL(EB), INTENT(INOUT) :: ZZ_GET(1:N_TRACKED_SPECIES)
REAL(EB) :: SUM_OTHER_SPECIES
INTEGER :: N_ZZ_MAX

! clip mass fractions
ZZ_GET=MAX(0._EB,MIN(1._EB,ZZ_GET))

! absorb all error in most abundant species
N_ZZ_MAX = MAXLOC(ZZ_GET,1)

SUM_OTHER_SPECIES = SUM(ZZ_GET) - ZZ_GET(N_ZZ_MAX)
ZZ_GET(N_ZZ_MAX) = 1._EB - SUM_OTHER_SPECIES

END SUBROUTINE GET_REALIZABLE_MF


SUBROUTINE GET_MASS_FRACTION(Z_IN,INDEX,Y_OUT)

! Y_OUT returns the mass fraction of species INDEX

INTEGER, INTENT(IN) :: INDEX
REAL(EB) :: Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: Y_OUT

Y_OUT = DOT_PRODUCT(Z2Y(INDEX,1:N_TRACKED_SPECIES),Z_IN)
Y_OUT = MIN(1._EB,MAX(0._EB,Y_OUT))

END SUBROUTINE GET_MASS_FRACTION


SUBROUTINE GET_MASS_FRACTION_ALL(Z_IN,Y_OUT)

! Y_OUT returns the mass fraction of all primitive species

REAL(EB) :: Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: Y_OUT(1:N_SPECIES)
INTEGER :: I

DO I=1,N_SPECIES
   Y_OUT(I) = DOT_PRODUCT(Z2Y(I,1:N_TRACKED_SPECIES),Z_IN)
ENDDO

Y_OUT = MIN(1._EB,MAX(0._EB,Y_OUT))

END SUBROUTINE GET_MASS_FRACTION_ALL


SUBROUTINE GET_MOLECULAR_WEIGHT(Z_IN,MW_OUT)

REAL(EB), INTENT(IN)  :: Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: MW_OUT

MW_OUT =  1._EB/DOT_PRODUCT(MWR_Z,Z_IN)

END SUBROUTINE GET_MOLECULAR_WEIGHT


SUBROUTINE GET_SPECIFIC_GAS_CONSTANT(Z_IN,RSUM_OUT)

REAL(EB) :: Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: RSUM_OUT

RSUM_OUT =  R0 * DOT_PRODUCT(MWR_Z,Z_IN)

END SUBROUTINE GET_SPECIFIC_GAS_CONSTANT


SUBROUTINE GET_SPECIFIC_HEAT(Z_IN,CP_OUT,TMPG)

INTEGER :: ITMP
REAL(EB), INTENT(IN) :: TMPG
REAL(EB) :: Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: CP_OUT

ITMP = MIN(5000,NINT(TMPG))
CP_OUT  = DOT_PRODUCT(CP_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)

END SUBROUTINE GET_SPECIFIC_HEAT


SUBROUTINE GET_SENSIBLE_ENTHALPY(Z_IN,H_S_OUT,TMPG)

INTEGER :: ITMP
REAL(EB), INTENT(IN) :: TMPG
REAL(EB) :: Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: H_S_OUT

ITMP = MIN(4999,INT(TMPG))

H_S_OUT = DOT_PRODUCT(H_SENS_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)+(TMPG-REAL(ITMP,EB))* &
          DOT_PRODUCT(H_SENS_Z(ITMP+1,1:N_TRACKED_SPECIES)-H_SENS_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)

END SUBROUTINE GET_SENSIBLE_ENTHALPY


SUBROUTINE GET_GIBBS_FREE_ENERGY(G_F_OUT,NU_IN,TMPG)

INTEGER :: ITMP
REAL(EB), INTENT(IN) :: TMPG,NU_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: G_F_OUT

ITMP = MIN(5000,NINT(TMPG))
G_F_OUT = 1.E6_EB*DOT_PRODUCT(G_F_Z(ITMP,1:N_TRACKED_SPECIES),NU_IN) ! convert from kJ/mol to J/kmol

END SUBROUTINE GET_GIBBS_FREE_ENERGY


SUBROUTINE GET_AVERAGE_SPECIFIC_HEAT(Z_IN,CPBAR_OUT,TMPG)

INTEGER :: ITMP
REAL(EB), INTENT(IN) :: TMPG
REAL(EB) :: Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: CPBAR_OUT

ITMP = MIN(5000,NINT(TMPG))
CPBAR_OUT = DOT_PRODUCT(CPBAR_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)

END SUBROUTINE GET_AVERAGE_SPECIFIC_HEAT


!SUBROUTINE GET_AVERAGE_SENSIBLE_SPECIFIC_HEAT(Z_IN,CP_AVG_OUT,TMPG)

!INTEGER :: ITMP
!REAL(EB), INTENT(IN) :: TMPG
!REAL(EB) :: Z_IN(1:N_TRACKED_SPECIES)
!REAL(EB), INTENT(OUT) :: CP_AVG_OUT
!
!ITMP = MIN(5000,NINT(TMPG))
!CP_AVG_OUT = DOT_PRODUCT(CP_AVG_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)
!
!END SUBROUTINE GET_AVERAGE_SENSIBLE_SPECIFIC_HEAT


SUBROUTINE GET_SENSIBLE_ENTHALPY_Z(N,TMPG,H_S)

INTEGER :: ITMP
REAL(EB), INTENT(IN) :: TMPG
INTEGER, INTENT(IN) :: N
REAL(EB), INTENT(OUT) :: H_S

ITMP = MIN(4999,INT(TMPG))
H_S = H_SENS_Z(ITMP,N)+(TMPG-REAL(ITMP,EB))*(H_SENS_Z(ITMP+1,N)-H_SENS_Z(ITMP,N))

END SUBROUTINE GET_SENSIBLE_ENTHALPY_Z


SUBROUTINE GET_AVERAGE_SPECIFIC_HEAT_Z(N,TMPG,CPBAR_OUT)

INTEGER :: ITMP
REAL(EB), INTENT(IN) :: TMPG
INTEGER, INTENT(IN) :: N
REAL(EB), INTENT(OUT) :: CPBAR_OUT

ITMP = MIN(5000,NINT(TMPG))
CPBAR_OUT = CPBAR_Z(ITMP,N)

END SUBROUTINE GET_AVERAGE_SPECIFIC_HEAT_Z


SUBROUTINE GET_CONDUCTIVITY(Z_IN,K_OUT,TMPG)

REAL(EB), INTENT(IN) :: Z_IN(1:N_TRACKED_SPECIES), TMPG
REAL(EB), INTENT(OUT) :: K_OUT
INTEGER :: ITMP

ITMP = MIN(5000,NINT(TMPG))
K_OUT = DOT_PRODUCT(K_RSQMW_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)/DOT_PRODUCT(Z_IN,RSQ_MW_Z)

END SUBROUTINE GET_CONDUCTIVITY


SUBROUTINE GET_SOLID_CONDUCTIVITY(K_OUT,TMP_S,OPT_MATL_INDEX,OPT_SURF_INDEX,OPT_RHO_IN,OPT_I_IN,OPT_J_IN,OPT_K_IN)

USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
INTEGER, INTENT(IN), OPTIONAL :: OPT_MATL_INDEX,OPT_SURF_INDEX,OPT_I_IN,OPT_J_IN,OPT_K_IN
REAL(EB), INTENT(IN), OPTIONAL :: OPT_RHO_IN(1:N_MATL)
REAL(EB), INTENT(IN) :: TMP_S
REAL(EB), INTENT(OUT) :: K_OUT
INTEGER :: N,NR,I_LOC,J_LOC,K_LOC
REAL(EB) :: VOLSUM
TYPE(MATERIAL_TYPE), POINTER :: ML=>NULL()
TYPE(SURFACE_TYPE), POINTER :: SF=>NULL()

!! K_OUT = 0.2_EB; RETURN ! PMMA debug

K_OUT = 0._EB
IF (PRESENT(OPT_MATL_INDEX)) THEN
   ML => MATERIAL(OPT_MATL_INDEX)
   IF (ML%K_S>0._EB) THEN
      K_OUT = ML%K_S
   ELSE
      NR = -NINT(ML%K_S)
      K_OUT = EVALUATE_RAMP(TMP_S,0._EB,NR)
   ENDIF
ELSEIF (PRESENT(OPT_SURF_INDEX)) THEN
   ! See FDS Tech Guide, Eqs. (7.22)-(7.25)
   SF => SURFACE(OPT_SURF_INDEX)
   K_OUT = 0._EB
   VOLSUM = 0._EB
   DO N=1,SF%N_MATL
      IF (OPT_RHO_IN(N)<TWO_EPSILON_EB) CYCLE
      ML => MATERIAL(SF%MATL_INDEX(N))
      VOLSUM = VOLSUM + OPT_RHO_IN(N)/ML%RHO_S
      IF (ML%K_S>0._EB) THEN
         K_OUT = K_OUT + OPT_RHO_IN(N)*ML%K_S/ML%RHO_S
      ELSE
         NR = -NINT(ML%K_S)
         K_OUT = K_OUT + OPT_RHO_IN(N)*EVALUATE_RAMP(TMP_S,0._EB,NR)/ML%RHO_S
      ENDIF
   ENDDO
   IF (VOLSUM>TWO_EPSILON_EB) K_OUT = K_OUT/VOLSUM
ENDIF
IF (K_OUT<TWO_EPSILON_EB) THEN
   IF (PRESENT(OPT_I_IN)) I_LOC = OPT_I_IN
   IF (PRESENT(OPT_J_IN)) J_LOC = OPT_J_IN
   IF (PRESENT(OPT_K_IN)) K_LOC = OPT_K_IN
   ! ! debug
   ! WRITE(LU_ERR,*) 'GET_SOLID_CONDUCTIVITY K_OUT: ',K_OUT
   ! WRITE(LU_ERR,*) 'GET_SOLID_CONDUCTIVITY VOLSUM: ',VOLSUM
   ! IF (PRESENT(OPT_I_IN)) WRITE(LU_ERR,*) 'GET_SOLID_CONDUCTIVITY I: ',I_LOC
   ! IF (PRESENT(OPT_J_IN)) WRITE(LU_ERR,*) 'GET_SOLID_CONDUCTIVITY J: ',J_LOC
   ! IF (PRESENT(OPT_K_IN)) WRITE(LU_ERR,*) 'GET_SOLID_CONDUCTIVITY K: ',K_LOC
   K_OUT = 10000._EB
ENDIF

END SUBROUTINE GET_SOLID_CONDUCTIVITY


SUBROUTINE GET_SOLID_RHOCBAR(RHOCBAR_OUT,TMP_S,OPT_MATL_INDEX,OPT_SURF_INDEX,OPT_RHO_IN)

USE MATH_FUNCTIONS, ONLY: EVALUATE_RAMP
INTEGER, INTENT(IN), OPTIONAL :: OPT_MATL_INDEX,OPT_SURF_INDEX
REAL(EB), INTENT(IN), OPTIONAL :: OPT_RHO_IN(1:N_MATL)
REAL(EB), INTENT(IN) :: TMP_S
REAL(EB), INTENT(OUT) :: RHOCBAR_OUT
INTEGER :: N,NR
TYPE(MATERIAL_TYPE), POINTER :: ML=>NULL()
TYPE(SURFACE_TYPE), POINTER :: SF=>NULL()

!! RHOCBAR_OUT = 2420000._EB; RETURN ! PMMA debug

RHOCBAR_OUT = 0._EB
IF (PRESENT(OPT_MATL_INDEX)) THEN
   ML => MATERIAL(OPT_MATL_INDEX)
   IF (ML%C_S>0._EB) THEN
      RHOCBAR_OUT = ML%RHO_S*ML%C_S
   ELSE
      NR = -NINT(ML%C_S)
      RHOCBAR_OUT = ML%RHO_S*EVALUATE_RAMP(TMP_S,0._EB,NR)
   ENDIF
ELSEIF (PRESENT(OPT_SURF_INDEX)) THEN
   ! See FDS Tech Guide, Eq. (7.24)
   SF => SURFACE(OPT_SURF_INDEX)
   DO N=1,SF%N_MATL
      IF (OPT_RHO_IN(N)<TWO_EPSILON_EB) CYCLE
      ML => MATERIAL(SF%MATL_INDEX(N))
      IF (ML%C_S>0._EB) THEN
         RHOCBAR_OUT = RHOCBAR_OUT + OPT_RHO_IN(N)*ML%C_S
      ELSE
         NR = -NINT(ML%C_S)
         RHOCBAR_OUT = RHOCBAR_OUT + OPT_RHO_IN(N)*EVALUATE_RAMP(TMP_S,0._EB,NR)
      ENDIF
   ENDDO
ENDIF
IF (RHOCBAR_OUT<=TWO_EPSILON_EB) RHOCBAR_OUT = 0.001_EB

END SUBROUTINE GET_SOLID_RHOCBAR


SUBROUTINE GET_SOLID_ABSORPTION_COEFFICIENT(KAPPA_OUT,SURF_INDEX,RHO_IN)

INTEGER, INTENT(IN) :: SURF_INDEX
REAL(EB), INTENT(IN) :: RHO_IN(1:N_MATL)
REAL(EB), INTENT(OUT) :: KAPPA_OUT
INTEGER :: N
REAL(EB) :: VOLSUM
TYPE(MATERIAL_TYPE), POINTER :: ML=>NULL()
TYPE(SURFACE_TYPE), POINTER :: SF=>NULL()

KAPPA_OUT = 0._EB
VOLSUM = 0._EB
SF => SURFACE(SURF_INDEX)
DO N=1,SF%N_MATL
   IF (RHO_IN(N)<=TWO_EPSILON_EB) CYCLE
   ML => MATERIAL(SF%MATL_INDEX(N))
   VOLSUM = VOLSUM + RHO_IN(N)/ML%RHO_S
   KAPPA_OUT = KAPPA_OUT + RHO_IN(N)*ML%KAPPA_S/ML%RHO_S
ENDDO
IF (VOLSUM>0._EB) KAPPA_OUT = KAPPA_OUT/VOLSUM

END SUBROUTINE GET_SOLID_ABSORPTION_COEFFICIENT


SUBROUTINE GET_SOLID_REFRACTIVE_INDEX(INDEX_OUT,SURF_INDEX,RHO_IN)

INTEGER, INTENT(IN) :: SURF_INDEX
REAL(EB), INTENT(IN) :: RHO_IN(1:N_MATL)
REAL(EB), INTENT(OUT) :: INDEX_OUT
INTEGER :: N
REAL(EB) :: VOLSUM
TYPE(MATERIAL_TYPE), POINTER :: ML=>NULL()
TYPE(SURFACE_TYPE), POINTER :: SF=>NULL()

INDEX_OUT = 0._EB
VOLSUM = 0._EB
SF => SURFACE(SURF_INDEX)
DO N=1,SF%N_MATL
   IF (RHO_IN(N)<=TWO_EPSILON_EB) CYCLE
   ML => MATERIAL(SF%MATL_INDEX(N))
   VOLSUM = VOLSUM + RHO_IN(N)/ML%RHO_S
   INDEX_OUT = INDEX_OUT + RHO_IN(N)*ML%REFRACTIVE_INDEX/ML%RHO_S
ENDDO
IF (VOLSUM>0._EB) INDEX_OUT = INDEX_OUT/VOLSUM

END SUBROUTINE GET_SOLID_REFRACTIVE_INDEX


SUBROUTINE GET_VISCOSITY(Z_IN,MU_OUT,TMPG)

INTEGER :: ITMP
REAL(EB), INTENT(IN)  :: TMPG,Z_IN(1:N_TRACKED_SPECIES)
REAL(EB), INTENT(OUT) :: MU_OUT

ITMP = MIN(5000,NINT(TMPG))
MU_OUT = DOT_PRODUCT(MU_RSQMW_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)/DOT_PRODUCT(Z_IN,RSQ_MW_Z)

END SUBROUTINE GET_VISCOSITY


SUBROUTINE GET_ENTHALPY(Z_IN,H_OUT,TMPG)

INTEGER :: ITMP
REAL(EB), INTENT(IN) :: TMPG
REAL(EB) :: Z_IN(1:N_TRACKED_SPECIES),DTMP
REAL(EB), INTENT(OUT) :: H_OUT

IF (TMPG>=5000._EB) THEN
   H_OUT = DOT_PRODUCT(CPBAR_Z(5000,1:N_TRACKED_SPECIES),Z_IN)*TMPG
ELSE
   ITMP = INT(TMPG)
   DTMP = TMPG-REAL(ITMP)
   H_OUT = DOT_PRODUCT(CPBAR_Z(ITMP,1:N_TRACKED_SPECIES),Z_IN)
   H_OUT = H_OUT+DTMP*(DOT_PRODUCT(CPBAR_Z(ITMP+1,1:N_TRACKED_SPECIES),Z_IN)-H_OUT)
   H_OUT = H_OUT*TMPG
ENDIF

END SUBROUTINE GET_ENTHALPY


REAL(EB) FUNCTION DRAG(RE,DRAG_LAW)

! drag coefficient

INTEGER, INTENT(IN) :: DRAG_LAW
REAL(EB), INTENT(IN) :: RE

SELECT CASE(DRAG_LAW)

   ! see J.P. Holman 7th Ed. Fig. 6-10
   CASE(SPHERE_DRAG)
      IF (RE<=TWO_EPSILON_EB) THEN
         DRAG = 100._EB
      ELSEIF (RE<=1._EB) THEN
         DRAG = 24._EB/RE
      ELSEIF (RE<1000._EB) THEN
         !!DRAG = 24._EB*(1._EB+0.15_EB*RE**0.687_EB)/RE ! see Crowe, Sommerfeld, Tsuji, 1998, Eq. (4.51)
         DRAG = 24._EB*(0.85_EB+0.15_EB*RE**0.687_EB)/RE ! matches Stokes drag at RE=1 (RJM)
      ELSEIF (RE>=1000._EB) THEN
         DRAG = 0.44_EB
      ENDIF

   ! see J.P. Holman 7th Ed. Fig. 6-9
   CASE(CYLINDER_DRAG)
      IF (RE<=1._EB) THEN
         DRAG = 10._EB/(RE**0.8_EB)
      ELSEIF (RE>1._EB .AND. RE<1000._EB) THEN
         DRAG = 10._EB*(0.6_EB+0.4_EB*RE**0.8_EB)/RE
      ELSEIF (RE>=1000._EB) THEN
         DRAG = 1._EB
      ENDIF

   CASE(USER_DRAG)
      DRAG = 1._EB ! PC%DRAG_COEFFICIENT set elsewhere

   CASE DEFAULT
      DRAG = 0._EB

END SELECT

END FUNCTION DRAG


REAL(EB) FUNCTION SURFACE_DENSITY(NM,MODE,WALL_INDEX,LAGRANGIAN_PARTICLE_INDEX,MATL_INDEX)

! Compute the surface density of a wall cell.
! If MODE=0, return kg/m2.
! If MODE=1, return kg/m3.

INTEGER, INTENT(IN) :: NM,MODE
INTEGER, INTENT(IN), OPTIONAL :: WALL_INDEX,LAGRANGIAN_PARTICLE_INDEX,MATL_INDEX
INTEGER :: I_GRAD=0,NWP,II2,N
REAL(EB) :: WGT,R_S(0:NWP_MAX)
TYPE(MESH_TYPE), POINTER :: M=>NULL()
TYPE(ONE_D_M_AND_E_XFER_TYPE), POINTER :: ONE_D
TYPE(SURFACE_TYPE), POINTER :: SF

M => MESHES(NM)

IF (PRESENT(WALL_INDEX)) THEN
   ONE_D => M%WALL(WALL_INDEX)%ONE_D
   SF => SURFACE(M%WALL(WALL_INDEX)%SURF_INDEX)
ELSEIF (PRESENT(LAGRANGIAN_PARTICLE_INDEX)) THEN
   ONE_D => M%LAGRANGIAN_PARTICLE(LAGRANGIAN_PARTICLE_INDEX)%ONE_D
   SF => SURFACE(LAGRANGIAN_PARTICLE_CLASS(M%LAGRANGIAN_PARTICLE(LAGRANGIAN_PARTICLE_INDEX)%CLASS_INDEX)%SURF_INDEX)
ELSE
   SURFACE_DENSITY = 0
   RETURN
ENDIF

IF (.NOT.SF%THERMALLY_THICK) THEN

   SURFACE_DENSITY = 0._EB

ELSE

   SELECT CASE(SF%GEOMETRY)
      CASE(SURF_CARTESIAN)    ; I_GRAD = 1
      CASE(SURF_CYLINDRICAL)  ; I_GRAD = 2
      CASE(SURF_SPHERICAL)    ; I_GRAD = 3
   END SELECT

   NWP = SUM(ONE_D%N_LAYER_CELLS)
   R_S(0:NWP) = SF%INNER_RADIUS + ONE_D%X(NWP) - ONE_D%X(0:NWP)

   SURFACE_DENSITY = 0._EB
   DO II2=1,NWP
      IF (MODE==0) THEN
         WGT = (R_S(II2-1)**I_GRAD-R_S(II2)**I_GRAD)/(REAL(I_GRAD,EB)*(SF%INNER_RADIUS+SF%THICKNESS)**(I_GRAD-1))
      ELSE
         WGT = (R_S(II2-1)**I_GRAD-R_S(II2)**I_GRAD)/(SF%INNER_RADIUS+SF%THICKNESS)**I_GRAD
      ENDIF
      IF (PRESENT(MATL_INDEX)) THEN
         SURFACE_DENSITY = SURFACE_DENSITY + ONE_D%MATL_COMP(MATL_INDEX)%RHO(II2)*WGT
      ELSE
         DO N=1,SF%N_MATL
            SURFACE_DENSITY = SURFACE_DENSITY + ONE_D%MATL_COMP(N)%RHO(II2)*WGT
         ENDDO
      ENDIF
   ENDDO

ENDIF

END FUNCTION SURFACE_DENSITY


SUBROUTINE PARTICLE_SIZE_DISTRIBUTION(DM,RR,CNF,CVF,NPT,GAMMA,SIGMA,DISTRIBUTION)

! Compute particle Cumulative Number Fraction (CNF) and Cumulative Volume Fraction (CVF)

USE MATH_FUNCTIONS, ONLY: IERFC
CHARACTER(LABEL_LENGTH), INTENT(IN) :: DISTRIBUTION
REAL(EB), INTENT(IN) :: DM,GAMMA,SIGMA
INTEGER, INTENT(IN) :: NPT
REAL(EB) :: SUM1,SUM2,DD1,DI,ETRM,GFAC,SFAC,DMIN,X1
INTEGER  :: J
REAL(EB), INTENT(OUT) :: RR(0:NPT),CNF(0:NPT),CVF(0:NPT)

RR(0)  = 0._EB
CNF(0) = 0._EB
SUM1   = 0._EB
SUM2   = 0._EB

X1     = IERFC(2._EB*CNF_CUTOFF)
DMIN   = MAX(DM*EXP(-X1*SQRT(2._EB)*SIGMA),0._EB)
DD1    = (-LOG(CNF_CUTOFF)/LOG(2._EB))**(1._EB/GAMMA)*DM
DD1    = (DD1-DMIN)/REAL(NPT,EB)
GFAC   = LOG(2._EB)*GAMMA*DD1/(DM**GAMMA)
SFAC   = DD1/(SQRT(TWOPI)*SIGMA)

INTLOOP: DO J=1,NPT
   DI = DMIN + (J-0.5_EB)*DD1
   RR(J) = 0.5_EB*DI
   IF ((DI<=DM .OR. DISTRIBUTION=='LOGNORMAL') .AND. DISTRIBUTION/='ROSIN-RAMMLER') THEN
      ETRM = EXP(-(LOG(DI/DM))**2/(2._EB*SIGMA**2))
      SUM1 = SUM1 + (SFAC/DI**4)*ETRM
      SUM2 = SUM2 + (SFAC/DI)*ETRM
   ELSE
      ETRM = EXP(-LOG(2._EB)*(DI/DM)**GAMMA)
      SUM1 = SUM1 + GFAC*DI**(GAMMA-4._EB)*ETRM
      SUM2 = 1._EB - ETRM
   ENDIF
   CNF(J) = SUM1
   CVF(J) = SUM2
ENDDO INTLOOP

CNF = CNF/SUM1
CVF = CVF/SUM2

END SUBROUTINE PARTICLE_SIZE_DISTRIBUTION


SUBROUTINE SPRAY_ANGLE_DISTRIBUTION(LON,LAT,LON_CDF,LAT_CDF,BETA,MU,SPRAY_ANGLE,DISTRIBUTION_TYPE,NPT)

INTEGER,INTENT(IN) :: NPT
REAL(EB),INTENT(OUT) :: LON_CDF(0:NPT),LON(0:NPT),LAT(0:NPT),LAT_CDF(0:NPT,0:NPT)
REAL(EB),INTENT(IN) :: BETA,MU,SPRAY_ANGLE(2,2)
CHARACTER(LABEL_LENGTH),INTENT(IN) :: DISTRIBUTION_TYPE
INTEGER :: I,J
REAL(EB) :: dLON,dLAT,PDF(0:NPT,0:NPT),THETA_MAX,THETA_MIN

THETA_MAX=MAXVAL(SPRAY_ANGLE)
THETA_MIN=MINVAL(SPRAY_ANGLE)

dLAT=(THETA_MAX-THETA_MIN)/REAL(NPT)
dLON=2._EB*PI/REAL(NPT)
!Discretize latitude and longtitude
LAT=(/ (THETA_MIN+I*dLAT,I=0,NPT) /)
LON=(/ (0._EB+I*dLON, I=0,NPT) /)

! SPRAY_ANGLE May be different in X and Y directions
! i.e spray angle is dependent on longtitude
! SPRAY_AGLE_MIN and MAX form ellipses with semi-axes defined by
! SPRAY_ANGLE(1,1:2) and SPRAY_ANGLE(2,1:2) respectively

DO I=0,NPT
  THETA_MIN=SPRAY_ANGLE(1,1)*SPRAY_ANGLE(1,2)
  IF(THETA_MIN>0._EB) THEN
     THETA_MIN=THETA_MIN/SQRT((SPRAY_ANGLE(1,2)*COS(LON(I)))**2+(SPRAY_ANGLE(1,1)*SIN(LON(I)))**2)
  ENDIF
  THETA_MAX=SPRAY_ANGLE(2,1)*SPRAY_ANGLE(2,2)
  THETA_MAX=THETA_MAX/SQRT((SPRAY_ANGLE(2,2)*COS(LON(I)))**2+(SPRAY_ANGLE(2,1)*SIN(LON(I)))**2)
  SELECT CASE(DISTRIBUTION_TYPE)
   CASE("TRIANGLE")
      DO J=0,NPT
         IF(LAT(J)<THETA_MIN .OR. LAT(J)>THETA_MAX) THEN
           PDF(J,I)=0._EB
         ELSE
           IF(LON(I)<MU) THEN
              PDF(J,I)=2*(LAT(J)-THETA_MIN)/((THETA_MAX-THETA_MIN)*(THETA_MAX-MU))
           ELSE
              PDF(J,I)=2*(THETA_MAX-LAT(J))/((THETA_MAX-THETA_MIN)*(THETA_MAX-MU))
           ENDIF
         ENDIF
      ENDDO
   CASE("GAUSSIAN")
      DO J=0,NPT
        IF(LAT(J)<THETA_MIN .OR. LAT(J)>THETA_MAX) THEN
           PDF(J,I)=0._EB
        ELSE
           PDF(J,I)=exp(-BETA*((LAT(J)-MU)/(THETA_MAX-THETA_MIN))**2)
        ENDIF
      ENDDO
   CASE DEFAULT ! "UNIFORM"
       DO J=0,NPT
        IF(LAT(J)<THETA_MIN .OR. LAT(J)>THETA_MAX) THEN
           PDF(J,I)=0._EB
        ELSE
           PDF(J,I)=1._EB
        ENDIF
      ENDDO
   END SELECT
ENDDO


!
DO I=0,NPT
PDF(I,:)=PDF(I,:)*SIN(LAT(I))
ENDDO

LAT_CDF=0._EB
! Latitude distribution conditional on Longtitude
DO I=1,NPT
   LAT_CDF(I,:)=LAT_CDF(I-1,:)+0.5*(PDF(I,:)+PDF(I-1,:))*dLAT
ENDDO

! Marginal longtitude distribution
LON_CDF=0._EB
DO I=1,NPT
   LON_CDF(I)=LON_CDF(I-1)+0.5*(LAT_CDF(NPT,I-1)+LAT_CDF(NPT,I))*dLON
ENDDO

! Normalize marginal longtitude distribution
LON_CDF=LON_CDF/LON_CDF(NPT)
!Normalize conditional latitude distributions
DO I=1,NPT
   LAT_CDF(:,I)=LAT_CDF(:,I)/LAT_CDF(NPT,I)
ENDDO

END SUBROUTINE SPRAY_ANGLE_DISTRIBUTION


REAL(EB) FUNCTION FED(Y_IN,RSUM,FED_ACTIVITY)

! Returns the integrand of FED (Fractional Effective Dose) calculation.

REAL(EB), INTENT(IN) :: Y_IN(1:N_TRACKED_SPECIES),RSUM
INTEGER, INTENT(IN) :: FED_ACTIVITY
INTEGER  :: N
REAL(EB) :: Y_MF_INT, TMP_1
REAL(EB), DIMENSION(3) :: CO_FED_FAC
!                at rest           light work(default) heavy work
DATA CO_FED_FAC /0.70486250E-5_EB, 2.7641667E-5_EB,    8.2925E-5_EB/

! All equations from D.A. Purser, Sec. 2, Chap. 6, SFPE Handbook, 4th Ed.
! Note: Purser uses minutes, here dt is in seconds. Conversion at the end of the function.
! Total FED dose:
! FED_dose = (FED_LCO + FED_LCN + FED_LNOx + FLD_irr)*FED_VCO2 + FED_LO2;

FED = 0._EB

! Carbon monoxide (CO)
!          at rest    light work heavy work
! RMV_FED /8.5_EB,    25.0_EB,   50.0_EB /
! D_FED   /40.0_EB,   30.0_EB,   20.0_EB /
! RMV/D   /0.2125_EB, 0.8333_EB, 2.5_EB/
!
! FED_LCO = (3.317E-5 * (C_CO)^1.036 * RMV * (dt/60)) / D;
!   with RMV=25 [l/min], D=30 [%] COHb concentration at incapacitation and C_CO in ppm
!
IF (CO_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,CO_INDEX,Y_MF_INT)
   TMP_1 = SPECIES(CO_INDEX)%RCON*Y_MF_INT*1.E6_EB/RSUM
   ! FED   = 2.764E-5_EB*TMP_1**(1.036_EB)
   FED   = CO_FED_FAC(FED_ACTIVITY)*TMP_1**(1.036_EB)
ENDIF

! Nitrogen oxides (NOx, here NO + NO2)
! FED_LNOx = C_NOx/1500 * (dt/60);
!   with C_NOx = C_NO + C_NO2, all in ppm
TMP_1 = 0._EB
IF (NO_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,NO_INDEX,Y_MF_INT)
   TMP_1 = SPECIES(NO_INDEX)%RCON*Y_MF_INT/RSUM
ENDIF
IF (NO2_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,NO2_INDEX,Y_MF_INT)
   TMP_1 = TMP_1 + SPECIES(NO2_INDEX)%RCON*Y_MF_INT/RSUM
ENDIF
IF (TMP_1 > 0._EB) FED = FED + TMP_1/0.001500_EB

! Cyanide
! FED_LCN = (exp(C_CN/43)/220 - 0.0045) * (dt/60);
!   with C_CN = C_HCN - C_NOx, all in ppm
IF (HCN_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,HCN_INDEX,Y_MF_INT)
   TMP_1 = SPECIES(HCN_INDEX)%RCON*Y_MF_INT/RSUM - TMP_1
   IF (TMP_1 > 0._EB) FED = FED + (Exp(TMP_1/0.000043_EB)/220.0_EB-0.00454545_EB)
ENDIF

! Irritants
! FLD_irr = (C_HCl/F_HCl + C_HBr/F_HBr + C_HF/F_HF + C_SO2/F_SO2 + C_NO2/F_NO2 + C_C3H4O/F_C3H4O + C_CH2O/F_CH2O) * (dt/60);
!   all in ppm
TMP_1 = 0._EB
DO N=1,N_SPECIES
   IF (SPECIES(N)%FLD_LETHAL_DOSE > 0._EB) THEN
      Call GET_MASS_FRACTION(Y_IN,N,Y_MF_INT)
      TMP_1 = TMP_1 + SPECIES(N)%RCON*Y_MF_INT/RSUM / SPECIES(N)%FLD_LETHAL_DOSE
   ENDIF
ENDDO
FED = FED + TMP_1

! Carbon dioxide (CO2) induced hyperventilation:
! FED_VCO2 = exp(0.1903*C_CO2/1E4 + 2.0004)/7.1;
!   C_CO2 in ppm
IF (CO2_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,CO2_INDEX,Y_MF_INT)
   TMP_1 = SPECIES(CO2_INDEX)%RCON*Y_MF_INT/RSUM
   If ( TMP_1 > 0.0_EB ) FED = FED * Exp( 0.1903_EB*TMP_1*100.0_EB + 2.0004_EB )/7.1_EB
ENDIF

! Low oxygen (O2)
! FED_LO2 = 1/exp(8.13 - 0.54*(0.209 - C_O2/1E6)) * (dt/60);
!   C_O2 in ppm
IF (O2_INDEX > 0) THEN
   Call GET_MASS_FRACTION(Y_IN,O2_INDEX,Y_MF_INT)
   TMP_1 = SPECIES(O2_INDEX)%RCON*Y_MF_INT/RSUM
   IF ( TMP_1 < 0.20_EB ) FED = FED + 1.0_EB  / Exp(8.13_EB-0.54_EB*(20.9_EB-100.0_EB*TMP_1))
ENDIF

! Convert the FED integrand for minutes.
FED = FED / 60._EB

END FUNCTION FED


REAL(EB) FUNCTION FIC(Y_IN,RSUM)
! Returns FIC (Fractional Incapacitating Concentration)

REAL(EB), INTENT(IN) :: Y_IN(1:N_TRACKED_SPECIES),RSUM
REAL(EB) :: Y_MF_INT
INTEGER  :: N

FIC = 0._EB
DO N=1,N_SPECIES
   IF (SPECIES(N)%FIC_CONCENTRATION > 0._EB) THEN
      Call GET_MASS_FRACTION(Y_IN,N,Y_MF_INT)
      FIC = FIC + SPECIES(N)%RCON*Y_MF_INT/RSUM / SPECIES(N)%FIC_CONCENTRATION
   ENDIF
ENDDO

END FUNCTION FIC


REAL(EB) FUNCTION WATER_VAPOR_MASS_FRACTION(HUMIDITY,TEMP,PZONE)

! Compute the water vapor mass fraction given the relative humidity and temperature

REAL(EB), INTENT(IN) :: TEMP,HUMIDITY,PZONE
REAL(EB) :: X_SAT,DHOR,P_RATIO
REAL(EB),PARAMETER :: T_BOIL=373.15_EB

DHOR = H_V_H2O(NINT(TEMP))*MW_H2O/R0
P_RATIO = P_STP/PZONE

X_SAT  = MIN(1._EB,P_STP/PZONE*EXP(DHOR*(1._EB/T_BOIL-1._EB/TEMP)))
WATER_VAPOR_MASS_FRACTION = HUMIDITY*0.01_EB*X_SAT/(MW_AIR/MW_H2O+(1._EB-MW_AIR/MW_H2O)*X_SAT)

END FUNCTION WATER_VAPOR_MASS_FRACTION


REAL(EB) FUNCTION RELATIVE_HUMIDITY(Y_H2O,TEMP,PZONE)

! Compute the relative humidity given the water vapor mass fraction and temperature

REAL (EB), INTENT(IN) :: TEMP,Y_H2O,PZONE
REAL (EB) :: X_SAT,X_H2O,DHOR
REAL(EB),PARAMETER :: T_BOIL=373.15_EB

DHOR = H_V_H2O(MIN(373,NINT(TEMP)))*MW_H2O/R0

X_SAT  = MIN(1._EB,P_STP/PZONE*EXP(DHOR*(1._EB/T_BOIL-1._EB/TEMP)))
X_H2O = Y_H2O*MW_AIR/(MW_H2O-Y_H2O*(MW_H2O-MW_AIR))
RELATIVE_HUMIDITY = 100._EB * X_H2O / X_SAT

END FUNCTION RELATIVE_HUMIDITY


REAL(EB) FUNCTION LES_FILTER_WIDTH_FUNCTION(DX,DY,DZ)
   REAL(EB), INTENT(IN):: DX,DY,DZ
   IF (TWO_D) THEN
      LES_FILTER_WIDTH_FUNCTION = SQRT(DX*DZ)
   ELSE
      LES_FILTER_WIDTH_FUNCTION = (DX*DY*DZ)**ONTH
   ENDIF
END FUNCTION LES_FILTER_WIDTH_FUNCTION


REAL(EB) FUNCTION GET_POTENTIAL_TEMPERATURE(TMP_IN,Z_IN)

USE MATH_FUNCTIONS, ONLY : EVALUATE_RAMP
REAL(EB), INTENT(IN) :: TMP_IN,Z_IN
REAL(EB) :: PP,DUMMY

PP = EVALUATE_RAMP(Z_IN,DUMMY,I_RAMP_P0_Z)
GET_POTENTIAL_TEMPERATURE = TMP_IN*(1.E5_EB/PP)**GM1OG  ! GM1OG = (GAMMA-1)/GAMMA = R/CP

END FUNCTION GET_POTENTIAL_TEMPERATURE


SUBROUTINE MONIN_OBUKHOV_SIMILARITY(Z,Z_0,L,U_STAR,THETA_STAR,THETA_0,U,TMP)

REAL(EB), INTENT(IN) :: Z,Z_0,L,U_STAR,THETA_STAR,THETA_0
REAL(EB), INTENT(OUT) :: U,TMP
REAL(EB), PARAMETER :: KAPPA=0.41_EB,P_0=100000._EB
REAL(EB) :: ZETA,PSI_M,PSI_H,THETA

IF (L>=0._EB) THEN
  PSI_M = -5._EB*Z/L
  PSI_H = PSI_M
ELSE
  ZETA = (1._EB-16._EB*Z/L)**0.25_EB
  PSI_M = 2._EB*LOG(0.5_EB*(1._EB+ZETA)) + LOG(0.5_EB*(1._EB+ZETA**2)) - 2._EB*ATAN(ZETA) + 0.5_EB*PI
  PSI_H = 2._EB*LOG(0.5_EB*(1._EB+ZETA**2))
ENDIF
U = (U_STAR/KAPPA)*(LOG(Z/Z_0)-PSI_M)
THETA = THETA_0 + (THETA_STAR/KAPPA)*(LOG(Z/Z_0)-PSI_H)
TMP = THETA*(P_0/(P_0-RHOA*GRAV*Z))**(-0.286_EB)

END SUBROUTINE MONIN_OBUKHOV_SIMILARITY


SUBROUTINE TURBULENT_VISCOSITY_INTERP1D(NU_OUT,NU_2,NU_3,DX_1,DX_2,DX_3)

! quadratic interpolation of the turbulent viscosity (MU_SGS) to the first off-wall gas phase cell center
! using 2nd and 3rd gas values
! assumes NU(z=0)=0, allows for non-uniform grid
!
!        NU_OUT         NU_2                 NU_3
! /////|    o    |        o        |           o           |
!         DX_1          DX_2                 DX_3

REAL(EB), INTENT(IN) :: NU_2,NU_3,DX_1,DX_2,DX_3
REAL(EB), INTENT(OUT) :: NU_OUT
REAL(EB) :: B,C,L_1,L_2,L_3

L_1 = 0.5_EB*DX_1
L_2 = DX_1 + 0.5_EB*DX_2
L_3 = DX_1 + DX_2 + 0.5_EB*DX_3

C = ( NU_3 - NU_2*L_3/L_2 ) / ( L_3**2 - L_2*L_3 )
B = ( NU_2 - C*L_2**2 ) / L_2

NU_OUT = MAX(0._EB,B*L_1 + C*L_1**2)

END SUBROUTINE TURBULENT_VISCOSITY_INTERP1D


END MODULE PHYSICAL_FUNCTIONS


MODULE TRAN

! Coordinate transformation functions

USE PRECISION_PARAMETERS
IMPLICIT NONE
TYPE TRAN_TYPE
   REAL(EB), POINTER, DIMENSION(:,:) :: C1=>NULL(),C2=>NULL(),C3=>NULL(),CCSTORE=>NULL(),PCSTORE=>NULL()
   INTEGER, POINTER, DIMENSION(:,:) :: IDERIVSTORE=>NULL()
   INTEGER NOC(3),ITRAN(3),NOCMAX
END TYPE TRAN_TYPE
TYPE (TRAN_TYPE), ALLOCATABLE, TARGET, DIMENSION(:) :: TRANS


CONTAINS


REAL(EB) FUNCTION G(X,IC,NM)

! Coordinate transformation function

REAL(EB), INTENT(IN) :: X
INTEGER, INTENT(IN)  :: IC,NM
INTEGER :: I,II,N
TYPE (TRAN_TYPE), POINTER :: T=>NULL()

T => TRANS(NM)

N = T%NOC(IC)
IF (N==0) THEN
   G = X
   RETURN
ENDIF

SELECT CASE(T%ITRAN(IC))
   CASE(1)
      G = 0._EB
      DO I=1,N+1
         G = G + T%C1(I,IC)*X**I
      ENDDO
   CASE(2)
      ILOOP: DO I=1,N+1
         II = I
         IF (X<=T%C1(I,IC)) EXIT ILOOP
      ENDDO ILOOP
      G = T%C2(II-1,IC) + T%C3(II,IC)*(X-T%C1(II-1,IC))
END SELECT

END FUNCTION G


REAL(EB) FUNCTION GP(X,IC,NM)

! Derivative of the coordinate transformation function

REAL(EB), INTENT(IN) :: X
INTEGER, INTENT(IN)  :: IC,NM
INTEGER :: I,II,N
TYPE (TRAN_TYPE), POINTER :: T=>NULL()

T => TRANS(NM)
N =  T%NOC(IC)
IF (N==0) THEN
   GP = 1._EB
   RETURN
ENDIF

SELECT CASE(T%ITRAN(IC))
   CASE(1)
      GP = 0._EB
      DO I=1,N+1
         GP = GP + I*T%C1(I,IC)*X**(I-1)
      ENDDO
   CASE(2)
      ILOOP: DO I=1,N+1
         II = I
         IF (X<=T%C1(I,IC)) EXIT ILOOP
      ENDDO ILOOP
      GP = T%C3(II,IC)
   CASE DEFAULT
      GP = 0._EB
END SELECT

END FUNCTION GP


REAL(EB) FUNCTION GINV(Z,IC,NM)

! Inverse of the coordinate transformation function

REAL(EB) :: GF
INTEGER :: N,IT,II,I
REAL(EB), INTENT(IN) :: Z
INTEGER, INTENT(IN)  :: IC,NM
TYPE (TRAN_TYPE), POINTER :: T=>NULL()

T => TRANS(NM)
GINV = Z
N = T%NOC(IC)
IF (N==0) RETURN

SELECT CASE(T%ITRAN(IC))
   CASE(1)
      LOOP1: DO IT=1,10
         GF = G(GINV,IC,NM)-Z
         IF (ABS(GF)<0.00001_EB) EXIT LOOP1
         GINV = GINV - GF/GP(GINV,IC,NM)
      ENDDO LOOP1
   CASE(2)
      ILOOP: DO I=1,N+1
         II = I
         IF (Z<=T%C2(I,IC)) EXIT ILOOP
      ENDDO ILOOP
      GINV = T%C1(II-1,IC) + (Z-T%C2(II-1,IC))/T%C3(II,IC)
END SELECT

END FUNCTION GINV


SUBROUTINE GET_IJK(X,Y,Z,NM,XI,YJ,ZK,I,J,K)

! Given the point (X,Y,Z) in mesh NM, GET_IJK returns the continuous and discrete form of the cell indices, (XI,YJ,ZK) and (I,J,K).
! For example, if the mesh is the unit cube with 10 cm cells, for the point (X,Y,Z)=(0.23,0.46,0.66), GET_IJK would return
! (XI,JK,ZK)=(2.3,4.6,6.6) and (I,J,K)=(3,5,7).

USE MESH_VARIABLES
REAL(EB),INTENT(IN) :: X,Y,Z
INTEGER,INTENT(IN) :: NM
REAL(EB), INTENT(OUT) :: XI,YJ,ZK
INTEGER, INTENT(OUT) :: I,J,K

XI = MESHES(NM)%CELLSI(FLOOR((X-MESHES(NM)%XS)*MESHES(NM)%RDXINT))
YJ = MESHES(NM)%CELLSJ(FLOOR((Y-MESHES(NM)%YS)*MESHES(NM)%RDYINT))
ZK = MESHES(NM)%CELLSK(FLOOR((Z-MESHES(NM)%ZS)*MESHES(NM)%RDZINT))
I = FLOOR(XI+1._EB)
J = FLOOR(YJ+1._EB)
K = FLOOR(ZK+1._EB)

END SUBROUTINE GET_IJK


END MODULE TRAN


MODULE OPENMP

! Module for OpenMP check

USE GLOBAL_CONSTANTS, ONLY : OPENMP_AVAILABLE_THREADS, OPENMP_USED_THREADS, OPENMP_USER_SET_THREADS, USE_OPENMP
!$ USE OMP_LIB
IMPLICIT NONE
PUBLIC OPENMP_INIT, OPENMP_SET_THREADS, OPENMP_PRINT_STATUS

CONTAINS

! set the control flag USE_OPENMP if OpenMP is used
SUBROUTINE OPENMP_INIT

!$OMP PARALLEL
!$OMP MASTER
!$ IF (OMP_GET_NUM_THREADS() /= 0) THEN
!$    USE_OPENMP = .TRUE.
!$    OPENMP_AVAILABLE_THREADS = OMP_GET_NUM_THREADS()
!$ ENDIF
!$OMP END MASTER
!$OMP BARRIER
!$OMP END PARALLEL

END SUBROUTINE OPENMP_INIT

! change the number of OpenMP threads if set by the user in the input file
SUBROUTINE OPENMP_SET_THREADS

!$IF (OPENMP_USER_SET_THREADS .EQV. .TRUE.) THEN
!$  IF (OPENMP_AVAILABLE_THREADS .NE. OPENMP_USED_THREADS) THEN
!$      CALL OMP_SET_NUM_THREADS(OPENMP_USED_THREADS)
!$  END IF
!$ELSE
!$  OPENMP_USED_THREADS = OPENMP_AVAILABLE_THREADS
!$  CALL OMP_SET_NUM_THREADS(OPENMP_USED_THREADS)
!$END IF

END SUBROUTINE OPENMP_SET_THREADS

! print OpenMP status
SUBROUTINE OPENMP_PRINT_STATUS
  USE GLOBAL_CONSTANTS, ONLY : LU_ERR, MYID, N_MPI_PROCESSES, VERBOSE
  INTEGER :: THREAD_ID

  !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(THREAD_ID)

  THREAD_ID = 0
  !$ THREAD_ID = OMP_GET_THREAD_NUM()

  !$OMP CRITICAL
  IF (USE_OPENMP .AND. OPENMP_USED_THREADS>1 .AND. VERBOSE) WRITE(LU_ERR,91) " OpenMP thread ",THREAD_ID," of ",&
     OPENMP_USED_THREADS-1," assigned to MPI process ",MYID," of ",N_MPI_PROCESSES-1
  IF (.NOT.USE_OPENMP .AND. VERBOSE) WRITE(LU_ERR,92) " MPI process ",MYID," of ",N_MPI_PROCESSES-1
  !$OMP END CRITICAL

  !$OMP END PARALLEL

  91 FORMAT(A,I3,A,I3,A,I6,A,I6)
  92 FORMAT(A,I6,A,I6)

END SUBROUTINE OPENMP_PRINT_STATUS

END MODULE OPENMP



MODULE MISC_FUNCTIONS

USE PRECISION_PARAMETERS
IMPLICIT NONE

CONTAINS


SUBROUTINE WRITE_SUMMARY_INFO(LU)

USE GLOBAL_CONSTANTS
USE COMP_FUNCTIONS, ONLY : GET_DATE
USE MPI
CHARACTER(LABEL_LENGTH) :: DATE
INTEGER, INTENT(IN) :: LU
INTEGER :: MPIVERSION,MPISUBVERSION,MPILIBLENGTH,IERR
CHARACTER(LEN=MPI_MAX_LIBRARY_VERSION_STRING) :: MPILIBVERSION

CALL GET_DATE(DATE)

WRITE(LU,'(/A/)')      ' Fire Dynamics Simulator'
WRITE(LU,'(A,A)')      ' Current Date     : ',TRIM(DATE)
WRITE(LU,'(A,A)')      ' Revision         : ',TRIM(GITHASH_PP)
WRITE(LU,'(A,A)')      ' Revision Date    : ',TRIM(GITDATE_PP)
#ifdef COMPVER_PP
WRITE(LU,'(A,A)')      ' Compiler         : ',TRIM(COMPVER_PP)
#endif
WRITE(LU,'(A,A/)')     ' Compilation Date : ',TRIM(BUILDDATE_PP)
                      WRITE(LU,'(A,I6)')  ' MPI Enabled;    Number of MPI Processes:  ',N_MPI_PROCESSES
IF (.NOT. USE_OPENMP) WRITE(LU,'(A)')     ' OpenMP Disabled'
IF (USE_OPENMP)       WRITE(LU,'(A,I6)')  ' OpenMP Enabled; Number of OpenMP Threads: ',OPENMP_AVAILABLE_THREADS

CALL MPI_GET_LIBRARY_VERSION(MPILIBVERSION,MPILIBLENGTH,IERR)
CALL MPI_GET_VERSION(MPIVERSION,MPISUBVERSION,IERR)
WRITE(LU,'(/A,I1,A,I1)') ' MPI version: ',MPIVERSION,'.',MPISUBVERSION
WRITE(LU,'(A,A)') ' MPI library version: ',TRIM(MPILIBVERSION)

END SUBROUTINE WRITE_SUMMARY_INFO


SUBROUTINE SEARCH_CONTROLLER(NAME,CTRL_ID,DEVC_ID,DEVICE_INDEX,CONTROL_INDEX,INPUT_INDEX)

USE DEVICE_VARIABLES, ONLY: DEVICE,N_DEVC
USE CONTROL_VARIABLES, ONLY: CONTROL,N_CTRL
USE COMP_FUNCTIONS, ONLY: SHUTDOWN
CHARACTER(MESSAGE_LENGTH) :: MESSAGE
CHARACTER(*), INTENT (IN) :: NAME,CTRL_ID,DEVC_ID
INTEGER :: I, DEVICE_INDEX,CONTROL_INDEX
INTEGER , INTENT(IN) :: INPUT_INDEX

! There cannot be both a device and controller for any given entity

IF (DEVC_ID /= 'null' .AND. CTRL_ID /='null') THEN
   WRITE(MESSAGE,'(A,A,1X,I3,A)')  'ERROR: ',TRIM(NAME),INPUT_INDEX,' has both a device (DEVC) and a control (CTRL) specified'
   CALL SHUTDOWN(MESSAGE,PROCESS_0_ONLY=.FALSE.)
ENDIF

! Search for device

IF (DEVC_ID /= 'null') THEN
   DO I=1,N_DEVC
      IF (DEVICE(I)%ID==DEVC_ID) THEN
         DEVICE_INDEX = I
         RETURN
      ENDIF
   ENDDO
   WRITE(MESSAGE,'(A,A,A)')  'ERROR: DEVICE ',TRIM(DEVC_ID),' does not exist'
   CALL SHUTDOWN(MESSAGE)
ENDIF

! Search for controller

IF (CTRL_ID /= 'null') THEN
   DO I=1,N_CTRL
      IF (CONTROL(I)%ID==CTRL_ID) THEN
         CONTROL_INDEX = I
         RETURN
      ENDIF
   ENDDO
   WRITE(MESSAGE,'(A,A,A)')  'ERROR: CONTROL ',TRIM(CTRL_ID),' does not exist'
   CALL SHUTDOWN(MESSAGE)
ENDIF

END SUBROUTINE SEARCH_CONTROLLER


END MODULE MISC_FUNCTIONS
